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i-^^Disorder and long-range interactions are two of the key components that make material failure an 
^interesting playfield for the application of statistical mechanics. The cornerstone in this respect has 
een lattice models of the fracture in which a network of elastic beams, bonds or electrical fuses 
Cwith random failure thresholds are subject to an increasing external load. These models describe 
|ion a qualitative level the failure processes of real, brittle or quasi-brittle materials. This has been 
4-Jparticularly important in solving the classical engineering problems of material strength: the size 
C^ependence of maximum stress and its sample to sample statistical fluctuations. At the same time, 
"^J^attice models pose many new fundamental questions in statistical physics, such as the relation be- 
tween fracture and phase transitions. Experimental results point out to the existence of an intriguing 
■^^rackling noise in the acoustic emission and of self-affine fractals in the crack surface morphology. 
^5^ecent advances in computer power have enabled considerable progress in the understanding of such 
rHmodels. Among these partly still controversial issues, are the scaling and size effects in material 
I strength and accumulated damage, the statistics of avalanches or bursts of microfailures, and the 
'"^morphology of the crack surface. Here we present an overview of the results obtained with lattice 
^Jnodels for fracture, highlighting the relations with statistical physics theories and more conventional 
^fracture mechanics approaches. 
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1 Introduction 

Understanding how materials fracture is a fundamental problem of science and 
engineering even today, although the effects of geometry, loading conditions, 
and material characteristics on the strength of materials have been empirically 
investigated from antique times. Probably the first systematic mathematical 
analysis of the problem dates back to the pioneering works of Leonardo da 
Vinci [1] and Galileo Galilei [2]. In Leonardo's notebooks one finds an inter- 
esting description of a tension test on metal wires [1]. Attaching a slowly filling 
sand bag to an end of the wire, Leonardo noticed that longer wires failed ear- 
lier (see Fig. A simple analysis based on continuum mechanics, however, 
would suggest that wires of the same section should carry the same stress and 
hence display equal strength. This lead some later commentators to hypothe- 
size that when Leonardo wrote smaller length, what he really meant was larger 
diameter [3]. As noted in Ref. [4], there is a simpler explanation: continuum 
mechanics treats the material as homogeneous, but this is often very far from 
the truth, especially for metal wires forged in the Renaissance. Disorder has 
profound effects on material strength, since fracture typically nucleates from 
weaker spots like preexisting microcracks or voids. The longer the wire, the 
easier it is to find a weak spot, hence the strength will decrease with length. 



Figure 1. The setup devised by Leonardo da Vinci to test the strength of metal wires. One end of 
the wire is suspended and the other end is attached to a basket that was slowly filled with sand. 
When the wire breaks, one can weigh the sand collected in the basket, determining the strength. 

From the original notebook [1]. 

What we have just described is probably the first experimental verification 
of size effect in fracture. Size effects are rooted in disorder and hence neces- 
sitates a statistical treatment. In particular, different samples of the same 
material could generally be expected to display fluctuations in the strength. 
The importance of defects and fluctuations in fracture was first realized in 
the early twenties by the work of Griffith [5] and later formalized from the 
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statistical point of view by Wcibull [6]. An excellent discussion of the history 
of strength up to the middle of the twentieth century can be found in the 
book by Timoshenko [7] . At present the importance of a statistical treatment 
in determining the strength of materials is certainly widely accepted, but de- 
spite this fact standard textbooks on fracture mechanics devote little attention 
to this problem. In addition experimental evaluations of strength distribution 
and other statistical properties are relatively scarce in the literature. 

Besides the engineering aspects of the problem, in the last twenty years the 
statistical properties of fracture have attracted a wide interest in the statis- 
tical physics community [8-10]. In this context, fracture can be seen as the 
outcome of the irreversible dynamics of a long-range interacting, disordered 
system. Several experimental observations have revealed that fracture is a com- 
plex phenomenon, described by scale invariant laws [11]. Examples notably 
include the acoustic emission measured prior to fracture and the roughness 
of the fracture surface. The observation of scaling makes it attractive to dis- 
cuss breakdown processes in terms of statistical mechanics and in particular 
to search for analogies with phase transitions and critical phenomena. These 
give rise to probability distributions that bear resemblance with the empirical 
observations on fracture, in particular "power-law" ones. 

The statistical physics approach to fracture is often based on numerical 
simulation of lattice models, which allow for a relatively simple description 
of disorder and elasticity. The elastic medium is represented by a network of 
springs or beams. The disorder is modeled for instance by imposing random 
failure thresholds on the springs or by removing a fraction of the links. In this 
way the model retains the long range nature and the tcnsorial structure of the 
interactions. The local displacements can then be found by standard method 
for solving coupled linear equations. The lattice is loaded imposing appropriate 
boundary conditions and the fracture process can be followed step by step, in 
a series of quasi-cquilibria. The cornerstone in this respect has been for the 
last twenty years the Random Fuse Network (or Model, RFM) [12], a lattice 
model of the fracture of solid materials in which as a further key simplification 
vectorial elasticity has been substituted with a scalar field. 

As damage accumulates, one observes intriguing history effects: in brittle 
failure, microcracks grow via the failure of the weakest elements [6]. In ductile 
fracture, dislocations or plastic deformation accumulates, and in viscoelastic 
materials there is a variety of complex relaxation mechanisms of stress in ad- 
dition to irreversible failure. Meanwhile, the stress field evolves and develops a 
more and more inhomogeneous structure. This arises in the simplest case from 
microcrack interactions that can shield each other or act jointly to enhance 
the local stress at various points in the material. 

In addition to numerical simulations, lattice models are useful to obtain 
new theoretical ideas and approaches to the problem. The aim is to develop a 
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probabilistic description of fracture involving the interplay between "frozen" 
properties, such as microcracks, and "thermal" noise. The former is indepen- 
dent of time, quenched into a sample and is only affected by whatever external 
influences the sample is subjected to and their internal consequences (consider 
growing microcracks). The thermal noise comes from environmental effects and 
the coarse-grained fluctuations at the atomistic level. In this respect, so-called 
fiber bundle models represent a useful idealization and the basis for more com- 
plex theoretical developments [13,14]. The exactly solvable global load sharing 
model can be seen as a mean-field version of the RFM. As it is customary with 
phase transitions and critical phenomena, mean-field theory provides a quali- 
tative picture of the transition. In general, however, for elastic, brittle failure 
the critical stress is ruled by the behavior of the extremal values of the stress 
distribution, the highest stress(es) present in any sample. This means that the 
highest moments of the distribution are crucial, contrary to what happens for 
other macroscopic properties such as the elastic moduli, which depend only 
on average mean-field quantities. One should be aware of this problem when 
trying to extend the results of fiber bundle models to more realistic cases. 

At a more application-oriented level, network models provide a valid alterna- 
tive to standard finite element calculations, especially in cases where structural 
disorder and inhomogeneities are particularly strong [15]. For instance, models 
of this kind are extensively used to model concrete [16-20]. This approach can 
be relevant for the predictability of structural failure, the design of stronger 
materials and optimization of structural composition. 

The scope of this review is to present an introduction to statistical models 
of fracture. At the same time, the goal is to provide an overview of the current 
understanding gained with this approach. After at least twenty years of intense 
research in the subject, we feel the need for a summarizing work, complement- 
ing the early insight gained in the past with more reliable results obtained 
by large scale simulations, which became possible due to recent advances in 
computer power and novel computational algorithmic techniques. This allows 
to clarify some controversial issues, such as the scaling in the damage dynam- 
ics, the avalanches or bursts of micr of allures, and the scaling of the fracture 
surface in lattice models. The access to better data and experimental progress 
has also created new unsolved problems. This work gives us thus an opportu- 
nity to outline both experimental and theoretical prospects. A research path 
that will not be pursued here is that of dielectric breakdown [21], which has 
many parallels with the RFM and the role of disorder is likewise important 
(see e.g. [9,22]). 

The review is organized as follows: a brief introduction to elasticity and 
fracture mechanics is reported in Sec. [21 While this does not represent by any 
means an exhaustive account of this very extended field, it has the primary 
purpose to set the formalism and the basic notation to be used in the fol- 
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lowing. In SecOlwe discuss the experimental background that motivates the 
introduction of statistical models of fracture. Thus particular attention will 
be devoted to those aspects that call for a statistical approach, such as the 
strength distributions and the related size effects, the crack roughness, and the 
acoustic emission. In Sec. 0] we proceed with a detailed description of network 
models for fracture, starting from the scalar RFM up to more complex beam 
models. We will highlight the connections with standard finite elements meth- 
ods and atomistic calculations and discuss dynamic effects. This complements 
the main emphasis here, of general theories and scaling laws by pointing out 
what issues are left out on purpose from the models at hand. 

The main theoretical approaches are reviewed in Sec. El starting from fiber 
bundle models as an illustration. We next discuss the general properties of 
phase transitions and point out the connections for fracture. Finally, we discuss 
the role of disorder in fracture model and introduce percolation concepts. 
Sect.Elis devoted to a discussion of recent simulation results. We report results 
for the constitutive behavior, fracture strength and damage distributions, size 
effects, damage localization, crack roughness and avalanches. The discussion 
section (Sec. E| outlines current challenges and possible future trends. An 
Appendix includes a description of the main numerical algorithms. 



2 Elements of fracture mechanics 

Fracture mechanics describes the formation and propagation of cracks in terms 
of macroscopic field equations, starting from the theory of linear elasticity. It is 
not our aim to review here the broad field of linear elastic fracture mechanics 
(LEFM) . We provide only a short summary of the main concepts defining the 
basic terminology to be used in the following sections. In addition, we present 
a general discussion of the role of disorder in fracture. 



2.1 Theory of linear elasticity 

The theory of linear elasticity describes the small deformations of solid bodies 
under the action of external forces. Consider a solid at rest and parameterize 
it by a set of coordinates f. Assuming that there is a direct relation between 
the coordinates r' of the deformed solid with those of the undeformed solid, we 
define the displacement field as n = r ' — r . It is then convenient to introduce 
the the symmetric strain tensor as 




(1) 
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When a solid deforms under the action of external forces, interactions at 
the molecular scale provide restoring forces that tend to bring the solid back 
to its reference configuration. Each of the material particles of the solid is in 
equilibrium under the action of both internal force interactions and external 
applied forces. Since all the internal forces should cancel each other due the 
action-reaction principle, when the equilibrium of a volume element dV is 
considered, we are left only with surface tractions that act on the surface dS 
that contains the volume element dV and externally applied body forces (such 
as gravity) that act over the volume element dV. In other words, the only non- 
vanishing contribution to the force should come from the tractions applied at 
the surface dS and the external body forces that act over volume dV. Thus 
the external force vector F with components Fj is given by 



where bj are the components of body forces per unit volume, and tj are the 
components of the surface tractions (force per unit surface area). By nota- 
tion, the bold symbols refer to vector or tensor quantities. Based on Cauchy's 
hypothesis, the surface tractions are given by tj = njCij, where n is the 
unit outward normal to the surface dS and a is the Cauchy stress tensor 
with components aij. Using the divergence theorem, Eq. [2lcan be written as 
F = Jy[h + 'V ■ a] dV (where V is the gradient operator), and in component 
form as 



The equilibrium condition for the body implies F = 0, which in terms of the 
stress tensor can be written as 




(2) 






(4) 



When the body forces b are absent, Eq.0 reduces to 




(5) 



In the limit of small deformations, there is a linear relation between stress 
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and strain 

CTij = ''^Cijkieki, (6) 



where Cijki is the elastic moduh tensor. Eq. ^ is the most general linear ex- 
pression relating the strain and the stress tensor, but it becomes much simpler 
in the case of isotropic media: 

Cijki = ^Sij6ki + fJ.{Sik6ji + 6ii6jk) (7) 

where the parameters A and fi are known as Lame coefficients, and 6ij = 1 
when i = j and is otherwise. It is often convenient to write the strain tensor 
as the sum of a compressive strain, involving a volume change, and a shear 
strain, which does not imply any volume changes. Then the Hooke laws can 
be written as 

aij = K6ij ^ ekk + 2/i(eij - -6ij ^ Ckk)- (8) 



where K is the bulk modulus and is the shear modulus. Finally, another 
standard formulation of the Hooke law involves the Young modulus E = -^r^ 

and the Poisson ratio = ^ f ^^x+l^ 



= (] I ,.\ I + ^ ^.^ij 1^ ^kk I • (9) 

The equation for stress equilibrium for an isotropic elastic medium can be 
written, combining the Hooke law with Eq. ((21), in terms of the displacement 
field u the Lame coefficients A and ^ as 

(A + /x)V(V • n) + ^V^n = 0, (10) 

with appropriate boundary conditions. Note that we have used both () nota- 
tion and the bold face notation interchangeably to denote a vector quantity. 



2.2 Cracks in elastic media 



In mathematical terms, a crack can be seen as a boundary condition for the 
elastic strain field, with the additional complication that the precise shape of 
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Figure 2. An elliptic crack in a two dimensional medium. 

the crack is in principle unknown, and the shape depends on the stress field 
acting on it. To overcome this formidable problem, one can prescribe the crack 
geometry and then compute the stress field. This approximation can provide 
a useful guidance on the mechanism of crack propagation. 

A simple illustrative example of this procedure is provided by the solution of 
the stress field around an elliptic crack placed under tensile stress in an infinite 
two dimensional medium (see Fig. EJ- The solution of the elastic equation of 
equilibrium yields for instance the stress concentration at the edge of the ellipse 

K^ = ^ = l + ^. (11) 
a b 

where a is the nominal applied stress and amax is the maximum stress, and 
a and b are major and minor radii of the ellipse. In reality cracks are not 
elliptical but typically display a sharp tip. If we rewrite Eq. in terms of 
the curvature radius p = b'^/a, we see that Kt diverges as l/^/p as p ^ 0, and 
hence the concept of stress concentration becomes useless. 
It was first shown by Irwin that the stress around a crack tip follows 



a,j = ^^fi,(9,4>) (12) 



where Kt is the stress intensity factor and fij is a crack-tip function in terms 
of the angular variables 9 and 4>. The particular values of Kt and the func- 
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tional forms of fij depend on the particular foading conditions (i.e. the way 
stress is apphed to the body). In this respect, it is convenient to identify three 
fracture modes (see Fig. corresponding to tensile (I), shear (II) and tearing 
(III) conditions (i.e., Kt equals Kj in mode I, Kjj in mode II, and Km in 
mode III). A generic loading rule can typically be split into an appropriate 
combination of the three modes. The presence of singularity around the crack 
tip, however, is an artifact of linear elasticity. Real materials can not sustain 
arbitrarily large stresses and in general close to the crack tip deformation is 
ruled by non-linear and inelastic effects. In practice it is customary to define a 
fracture process zone (FPZ) as the region surrounding the crack where these 
effects take place, while far enough from the crack one recovers linear elastic 
behavior, that is the long-range decay of perturbed displacement fields. 

Knowing the stress field around a crack does not tell us directly how the 
fracture process would develop, but this information is necessary to build a 
model. A simple and clear argument to understand the stability of a crack 
was provided by Griffith almost a century ago [5]. The idea is to compute 
the energy needed to open a crack and see when this process becomes favor- 
able. Considering for simplicity a two dimensional geometry, the elastic energy 
change associated with a crack of length a is given by 

^ (13) 

where a is the applied stress. Forming a crack involves the rupture of atomic 
bonds, a process that has an energy cost proportional to the crack surface 

Esurf = ^a-f, (14) 

where 7 is the surface tension. Griffith noted that a crack will grow when this 
process leads to a decrease of the total energy E = E^^i + Eg^rf- Thus the 
criterion for crack growth can be expressed as 

dE 

_ = -vrcT^a/E^ + 27 < 0, (15) 
da 

which implies that under a stress a cracks of size a > Oc = 2'yE / {tto''^) are 
unstable. Equation (jlSf) can also be written as a^pKo, > \J2E^ and, noting 
that K = cjy^vra is the stress intensity factor of an infinite plate, we can 
reformulate the crack growth criterion in terms of a critical stress intensity 
factor Kc = \/2E'y, also known as the material toughness (i.e. the crack grows 
when K > Kc). 
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Mode I Mode II Mode HI 



Figure 3. The three modes of fracture. 

2.3 The role of disorder on material strength 

The influence of randomness on fracture can have many disguises. It is useful 
to leave aside inelastic effects for awhile, and consider only brittle or "quasi- 
brittle" processes. In this case, we can generalize the discussion of section 
using the concepts of long-range elastic fields, critical energy or crack size. 
There are several independent subcases that present different physics: one may 
consider a single crack and analyze the effect of disorder on material strength 
or, alternatively, one could study the dynamics of a whole crack population. 
In general, we can identify some key issues to address: 

(i) what is the effect of disorder on the stability of a single crack, and how 
does the Griffith argument get modified? 

(ii) in the presence of many cracks, how is the physics changed from the case 
of a single crack? 

(iii) what is the effect of disorder on the propagation of a single crack? 

(iv) what are the concepts, and statistical laws that operate in the case of many 
cracks of which no single one dominates over the others? 

These four fundamental questions also touch upon a number of engineer- 
ing mechanics and materials science issues, such as fracture toughness, crack 
arrest, defect populations, damage mechanics, and failure prediction. 

The first two issues from the above list deal with size scaling and the statis- 
tics of strength. Figure 0] illustrates the schematic advancement of a two- 
dimensional crack from a pre-existing flaw (or a notch in a test). The presence 
of variation in the material structure can be presented ideally - in brittle mate- 
rials - by locally varying elastic moduli, so that the stress-field becomes more 
complicated than the homogeneous LEFM solution. It can also be seen as a 
crack surface energy 7s(x) which has a fluctuating component. The issue of 
size effects amounts to computing the strength of a sample of linear size L. The 
answer of course depends on whether the sample has one dominating crack or 
more number of microcracks or flaws. 
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Figure 4. One-dimensional example of a crack propagating in a disordered environment. Included 
are crack arrest due to strong regions, deflection, branching and wandering. 

In the case of a single dominating defect, the Griffith argument can be eas- 
ily modified to account for a spatial variation of the elastic coefficients. The 
technical problem is to consider the distributions of the surface energy 7^ and 
elastic modulus or, even the presence of dislocations in the sample. As a 
result, one gets instead of the deterministic critical crack length (for constant 
E and 7s) a probability for the crack to be stable under a given external 
stress. In other words, we can characterize the problem by a strength distribu- 
tion which in the case of a single crack crucially depends on the microscopic 
randomness. For instance, Arndt and Nattermann [23] showed that an uncor- 
related random 7^ results in an exponential strength distribution. In general, 
however, 7^ could also display spatial correlations, as measured by the two 
point correlator (7s (2^)75 (y)), where the brackets denote a statistical average. 
Correlations in 7^ imply that an increased crack opening will cost energy that 
not only depends on 7^ at a single location but also depends on the presence 
of such correlations. 

2.4 Extreme statistics for independent cracks 

In the presence of many non-interacting defects, it becomes apparent that the 
physics is related to extreme statistics (see Fig.|2I). The simplest starting point 
is to find the largest defect - or the most dangerous one - in the system. Later 
on in the context of random fuse networks we will return to this argument in 
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more detail. If the critical crack fails instantaneously, there is no pre-critical 
crack growth, and the cracks do not interact, then this is a sufficient condition 
for determining the statistical strength properties. In other words, we want to 
consider a dilute gas of non-interacting defects, and find the worst one using 
a Lifshitz type argument [6]. 

Consider a collection of Nfi defects, and assign to each of these a failure 
strength o"j, where i = 1, . . . , Nd- The whole sample fails if one of these defects 
can not stand the applied stress. Thus for the cumulative strength distribution 
F{a, Nd) one has that 

1-Fia)=ll[l-P,{a)]. (16) 
1=1 

If the cracks are independent so that ai are random variables distributed 
according to Pi = P{o-), then it follows that 

1 - F{a) ~ ex.-p[-VpP{a)] (17) 

where p is the crack density and V the sample volume. The crucial ingredients 
here are the P{c7) and N^, since they rule both the asymptotic form of the 
distribution F as well as its form for finite N^. Of these two issues, more 
interesting is the behavior of F{a) in the large N^, limit. Then, the question 
boils down to the actual distribution of the P{cr). For narrow distributions, 
one has the Gumbel scaling form 

F{a) ~ 1 - exp [-V expp(c7)] (18) 

where the function p{a) follows from P{a) and depends on two things: the 
stress enhancement factor of a typical defect of linear size I, and the size 
distribution P{1) of the defects. For example, an exponential P(l) will thus 
result in this Gumbel scaling form. 

The other possibility is to use a distribution P{1) that makes the local 
stresses at crack tips to have a wide distribution as well. Using a length dis- 
tribution P(/) ~ /^^ changes the limiting distribution of the sample strength 
to a Weibull distribution. The reason is simple: assuming that the stress a at 
the crack tips grows as a power of crack length /, then the P{a) also becomes 
a power law distribution since P{1) is a power law distribution. As a result, 
the celebrated Weibull distribution [6] (or its 3-parameter version) arises: 

F{a) ~ 1 - exp i-Va'^/A'^) (19) 

where A is a normalization constant, and rn is the Weibull modulus, which 
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Figure 5. An example of a dilute defect population. These cracks each of some size I are assumed 
in the simplest case to have no influence on each other, and they have a linear size distribution P{1). 

can be related to the strength of a defect of size In particular, we have 
that if o" ~ then m = 2g. As usual, a small m implies a heterogeneous 

material/sample as g is then small as well. In any case, the distributions Eqs. 
(|18|) . (|19j) allow now to distinguish between microscopic details in experiments, 
and predict the size-scaling of average strength. The Gumbel-form is apt to 
give a logarithmic size-effect, but note that this is only true in the absence of 
notches or dominating defects (see Sec. I5.3.IT)) . 



2.5 Interacting cracks and damage mechanics 

The discourse on strength changes character if one allows for slow average crack 
growth. Here, we would like to focus on the case in which the fracture dynamics 
is continuous and slow, so that in spite of possibly very fast local developments, 
a sample can be considered in a quasi-equilibrium. This means that the sample 
internal state is either independent of time or evolves slowly with respect to 
the fast time-scale of the fracture process. Even in these conditions, treating 
explicitly the dynamics of interacting microcracks represents a formidable task: 
the growth of a microcrack can be either inhibited/screened or accelerated due 
to the presence of a neighboring microcrack. The outcome in general depends 
on the loading conditions, the crack geometry and other details, as discussed 
in the engineering mechanics literature. 

To avoid all these complications, one can disregard the dynamics of individ- 
ual microcracks and treat the damage at a coarse grained scale. This approach, 
known as "damage mechanics" [17,24,25], is valid when the interaction among 
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cracks is moderate. In the simplest case, one can consider the evolution of a 
scalar valued damage field D{x), describing the local variations of the Hooke 
law 



The idea behind damage mechanics is that all the complications of the fracture 
process can be encoded in the field D, that will then evolve according to a 
prescribed law, reflecting the internal microcracking in the material. 

Simple models of statistical fracture can be used to check some of the typical 
assumptions and possible outcomes in damage mechanics. This is in particular 
applicable in "scalar fracture", described by the electrical analogue of the 
RFM. In order to apply the coarse-graining procedure, implied by Eq. (I^U)) . 
one assumes the presence of a length-scale over which the damage self- averages, 
meaning that its fluctuations become much smaller than the average. This scale 
is also known as the Representative Volume Element (RVE) above which it is 
possible to treat the damage as a smoothly varying variable. 

Figure El depicts the typical scenario for damage accumulation. In this ex- 
ample we imagine a set of microcracks that will grow in a stable manner, and 
give rise to locally enhanced "damage" D[x). One may now naturally deflne 
a correlation length, by measuring e.g. the two-point correlation function of a 
scalar D via 



which should decay, for a fleld D which is uncorrelated over longer dis- 
tances, asymptotically exponentially and thus deflnes a correlation length ^ 
via Cd ~ exp(— r/^). An interesting issue is the "symmetry breaking" in 
damage accumulation: the direction of crack growth is mostly dictated by ge- 
ometry and loading conditions. Thus, as illustrated by Fig. El it is necessary to 
consider separately the directions perpendicular and parallel to the external 
stress, which would then imply two correlation lengths, and 

Given a fracture process, in which damage grows in a controlled fashion one 
can then - in principle, at least in simulations - follow the development of ^. 
This quantity will increase, making a sample effectively "smaller": as ^ (and 
the RVE size ) increases L/^ diminishes, where L is as usual the sample linear 
size. An additional complication comes from the damage anisotropy, since the 
two correlation lengths and will follow different dynamics. 

It is interesting to discuss what happens when one of the two becomes much 
larger than the other. In this case the sample becomes effectively much "larger" 
along the stress direction, becoming almost one-dimensional, which in other 
words indicates the presence of damage localization. The developments that 




(20) 



Cn{r) = {{{D{x + r) - {D)) - {D{x) - {D))f) 
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lead to and beyond this point are of course an interesting issue, with also 
materials science twists. Consider for instance the mechanics of fiber-reinforced 
composites: eventually one crack will dominate, but still the specimen may 
have a finite strength remaining thanks to fibers that bridge the crack and 
induce cohesion. 

In general, it is a complicated question as to what kind of evolution would 
rule {D), even if the possibly very important role of fluctuations is neglected. 
That is, as the ^ increases the damage inside the corresponding RVE should 
fluctuate more and more. Recall also the connection between D and the sample 
effective elastic modulus, which implies that the two follow similar dynamics. 
One can postulate several possible scenarios for the damage dynamics, in- 
cluding so-called "finite-time singularities" (see e.g. [26,27]). These imply a 
divergence of the strain under a linearly growing applied stress a (xt 



e~l/(te-t)^ (22) 

i.e. a power-law -like or pseudo-critical approach to the failure time tc and to 
the corresponding critical stress cTc oc tc- In the context of damage mechanics, 
Eq. (|22j) would result from a stress dependent damage law of the type 1 — D ~ 
(fjc — cr)"^. In cases with a dominant crack this law can be related to the crack 
growth rate as a function of a. Another possibility is that the damage displays 
an abrupt jump at ac- The two possibilities are depicted in Fig. El 

The statistical laws and geometry of damage accumulation can be studied 
both via laboratory experiments and the models that are the main scope of this 
review, and they often have also great importance for engineering applications. 
The typical engineering mechanics viewpoint of the previous discussion is to 
look at the rate of crack growth given just a single crack. The assumption is 
that the crack length a evolves according to 

ft = ^KTff (23) 

where -ft^e// is an effective stress intensity factor and the exponent m typically 
varies between 30 and 50 [28]. Eq. (|23|) . known as the "Charles law", is a 
way of generalizing the standard LEFM expression, in Eq. (|12j) . for which 
Ki,ff ~ a^/^. It implies that even a subcritical crack, in the sense of the 
Griffith argument, can slowly grow in some conditions such as creep or fatigue. 
Examples arise in non-brittle fracture, and in the presence of a diffuse crack 
population [29]. For i^e// ^ L one gets a power-law growth in time by simple 
integration [30] . The possible roles of disorder here is a question which has not 
been studied systematically. 
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Figure 6. An example of a "first-order" and a critical scenario for damage build-up, via 

Be// = SD(e) 



2.6 Fracture mechanics of rough cracks 

2.6.1 Crack dynamics in a disordered environment: self-affinity and 
anomalous scaling. While damage mechanics treats fracture in disordered 
media at a coarse grained scale, it is often necessary to treat the dynamics of 
a single macroscopic crack. A simple example is to consider a medium con- 
taining regions with sufficiently high 7^, locally, so as to stop crack growth. In 
these circumstances energetic considerations dictate that the crack path will 
deviate from the straight direction (see Fig. 0]). The simple case is the two- 
dimensional one, when the vertical coordinate of the crack tip, (where t 
is the position along the horizontal projection of the crack path) wanders with 
equal probability up and down. This process corresponds directly to a random 
walk or Brownian motion. This fact implies that the crack surface, correspond- 
ing to the trajectory of the crack tip, displays a characteristic roughness. In 
the following we consider the fracture surfaces /i(x, t) in various cases. The 
following sub-cases can be distinguished: 

(i) 2D crack: the crack line grows as a trail of a zero-dimensional point-like 
particle, the effective crack point. For example, consider the growth of a 
crack from a pre-made notch. 

(ii) 3D interfacial crack: a crack line propagates in a plane (xy), so that 
the z-direction can be neglected. If the average direction of propagation is 
in the y-direction, say, one is interested in the statistics of the rough line 
h{x, y) where y is the time-dependent average position, and h varies with 
X, in the y-direction. 

(iii) 3D crack: one has a crack line given by f[t). This leaves a trail as the time 
increases, which is a two-dimensional surface /i(x, y) where h is along the 
z-axis. In general, h is defined only for those (x, y) for which the projection 
of r on the xy-plane, {{x',y')), is such that that for x = x' , y' > y. One 
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may now study the projections of h to the (xy)- and ((x2;))-planes, defining 
thus the two components h\\{x,y) and h±{x,y) which he in the respective 
planes. 

The geometry of rough surfaces becomes particularly interesting when it 
turns out that they are self-affine fractals, or even more complicated versions 
thereof as discussed below. There are three core questions related to this kind 
of "roughening", in addition to the phenomenology of exponents and scaling 
relations, and the origins of such scaling: 

• What are the the smallest length- (and time-) scales on which one can see 
such behavior? 

• What is the largest scale, that of "saturation" , where the roughness becomes 
independent of scale? What leads to the saturation, in particular if one is 
concerned about the microstructure? 

• How does the roughness evolve? A crack is often induced in the laboratory 
from a notch, which is "flat" . The roughness of the crack will develop as new 
crack surface is created. Thus one must be careful in considering a posteriori 
surfaces from experiment or simulations. 

Let us now study the separate possibilities. Of these, case (ii) from above is 
conceptually the easiest one. Assume that the crack front advances on average 
with a constant velocity v, so that y = vt. The simplest quantity used to 
characterize the roughness is the global surface width 



which simply measures the standard deviation around the mean position. The 
global width is a quantity often considered in kinetic roughening of surfaces, for 
instance in molecular-beam epitaxy or in the propagation of slow combustion 
fronts in paper. In these cases, it happens that the dynamics of h{x, t) is 
ruled by a Langevin equation and one can expect the so-called Family- Vicsek- 
scaling [31]. This implies the presence of an unique spatial scale, the correlation 
length and an associated time-scale. In practice one writes 



Eq. (|25() describes the roughening of h from an initial condition, which is 
usually a flat one, although one can also consider the influence of an arbitrary 
profile h{x,t = 0). The initial transient is governed by the growth exponent 
P, until correlations build up, as ^ reaches L, and the roughness saturates at 
the value W ^ L^. This implies that the scaling function f{x) is a constant, 
except for x small. The self-affine scaling relates now /3, (, and the dynamic 




(24) 




(25) 
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exponent z through 

Hz = C. (26) 

Note that one can have further constraints, such as the so-called Galileian 
invariance for the Langevin equation, coupling the exponents (2 — = C in 
2D). 

If we consider separately the profiles /i||(x,y) and h±{x,y) (again, one may 
substitute y = vt), in the simplest imaginable case, the Family- Vicsek picture 
applies to both, separately, with six separate exponents zy, and P±, 
z±, C±, respectively, that are coupled by the relation Eq. [^Hl In section ESI 
we overview line models of cracks, according to scenario (iii), such that they 
might relate to this scenario. 

The existence of a correlation length begs the question of where it stems 
from in the crack growth dynamics, or in other words if one measures the ex- 
ponents experimentally what is the equation underlying the roughening pro- 
cess. This is apparent if we consider the generalization of the one-dimensional 
example of a Brownian "crack" trail to case (iii), which is the most relevant 
for experiments. The idea is to represent the opening of a crack by a moving 
one-dimensional line which leaves a two-dimensional surface h{x, y) trailing 
after itself. Thus in the intermediate fracture stage one has a line separating 
the untouched part and the crack - with two open faces. Considering the case 
slow fracture, one could then write equations of motion for the fluctuating 
components of the crack line, taking into account the pertinent physics as the 
presence of randomness, and the coupling due to elasticity of the fluctuations, 
both in-plane and out-of-plane. Such an approach also contains the idea that 
the physics in the bulk, as in the fracture process zone, can be neglected or 
incorporated again as a long-range elasticity along the crackline. Similar mod- 
els that include further complications such as dynamical overshoots along the 
crack front have also been studied to a great extent (see section . 

Scalings such as that implied by Eq. (|25|1 are to be understood in the av- 
erage sense. That is, one has to sample the stochastic process sufficiently in 
order to ensure the convergence of quantities. This is naturally a problem in 
fracture, where it may be difficult to collect enough samples. There are other 
measures that can be used to extract the exponents, such as the local width, 
the correlation functions, and the structure factor [32], but notice that the 
discussion assumes that the profile h is translation-invariant, which may be 
invalidated by the presence of boundary effects. 

An example thereof is the case (i) from above, which in the simplest imag- 
inable case might be equivalent to Brownian motion. Now, h(x) resembles the 
case of a propagating one-dimensional crack- front, but is actually different in 
several ways. First, the already- formed profile is frozen in time (unless one 
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has microscopic processes that change the free crack surfaces). Second, as in 
the Family- Vicsek picture of Eq. H25|) . one may assume that there are initial 
transients. Their character and the appropriate scaling picture would depend 
on the details. For instance, one can consider a Brownian motion confined to 
a maximum excursion, so that eventually the roughness would saturate. In 
reality, one should keep in mind that due to the presence of long range elastic 
fields, the dynamics of the crack tip is never a simple uncorrelated random 
walk in 2D. The interplay between elasticity and disorder gives in general rise 
to avalanches or stick-slip motion, as it is typical of non-linear driven systems. 
One may imagine that in 2D the crack is arrested at obstacles, slowing it down 
till a deviating path is found and the crack tip speeds up till the next obstacle. 
Such avalanches can be detected through acoustic emission and we will discuss 
this issue further in Sec. 13 

In real experiments, as in many that relate to case (iii) from above, one can 
also observe two complicated scenarios than Family- Vicsek scaling: anoma- 
lous scaling and multiscaling. To this end, we should consider more elaborate 
measures of rough surface profiles h{x), measured over a window of size I < L 
{x implies, that we will not distinguish the various crack- inspired scenarios). 
In full generality we can define the q-moments and assume that they scale as 
power-laws 



where 6h = h — h is the deviation of h inside a window of size I, and ^ is the 
correlation length at which Wq becomes a constant. The brackets imply again 
an ensemble average. The simplest case is encountered when all the g— order 
roughness exponents Cg are the same (i.e. Q = for all q), corresponding 
to self-afhne scaling. The opposite possibility is referred to as multiscaling. 
In the context of kinetic interface roughening multiscaling implies that addi- 
tional mechanisms create large height differences which are measured more 
sensitively by the higher q-moments (an example of this will be discussed in 



In terms of the local roughness, the Family- Vicsek scaling can be written as 



w, 



'g{l,t) = {{6h{xy))^/'' ~ l<- for / < C{t). 



(27) 



Sec.E23I). 




(28) 



with the scaling function W(x) ~ for x < 1 and a constant for x ^ 1. 
Notice that the exponent ( can also be obtained by the structure factor, or 
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spatial power spectrum 



S(k,t)^(|/l(k,t)|2) 




^-{2C+d) for J. ^ 
^(t)2C+'^ for A: < 1/^(0 



(29) 



where /i(k, t) denotes the spatial (d-dimensional) Fourier transform of h{x,t). 
Alternatively, one can use the height difference correlation function 



which for large separations scales as Ixp. 

If in Eq. H28|) the amplitude a{t) increases with time, anomalous scaling en- 
sues [33], and if C > 1 one talks about superroughness [33,34]. The quenched 
Edwards- Wilkinson equation in d = 1, describing the "depinning" of a line 
with local elasticity is a classical example of this [35-39]. Anomalous scal- 
ing usually implies that the roughness of a self-affine object (e.g. a crack) 
develops a power-law dependence on the system size L [33,34]. One should 
thus define two different exponents C and Cioc^ the first exponent describing 
the dependence of the global roughness on the system size, while the second 
encodes the scaling of the local roughness on the measuring window. In partic- 
ular, the local width scales as w{l) ~ /^'°'^L^~^'°% so that the global roughness 
w{L) scales as with ( > (\oc- Consequently, the power spectrum scales 
as S{k) ~ fc-(2C.oc+rf)^2(C-C,„c)_ Notice also that Eqs. EH and 03 reflect only 
the local roughness exponent (^loc. This implies that for fracture surfaces with 
two (independent) projections, one can define no less than four independent 
roughness exponents, and thus there could be eight exponents describing this 
case. 

An interesting issue here is how ordinary LEFM changes in the presence of 
a crack with a "rough" or self-affine geometry. The same problem can also 
be restated in the presence of general damage or a population of self-affine 
microcracks. All in all, there are at least three different possibilities. One can 
consider the stress-field in the proximity of a "fractal crack", but then one 
has to remember that Nature does not create a priori such objects - they 
grow because of some stable crack growth process. This takes place in the 
presence of elastic fields possibly modified from LEFM predictions due to the 
fractal geometry [40-42]. A particular version is to consider the Lame equation, 
Eq. H10|). in the presence of a given self-affine boundary. This can be solved in 
some cases to illuminate possible changes from LEFM, but is of course outside 
of fracture. Finally, a promising avenue is given by the work of Procaccia and 
collaborators [43-45], in which conformal mapping techniques developed in 
the context of Laplacian growth problems are combined with two-dimensional 




(30) 
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elasticity. This allows to grow relatively large cracks iteratively, accounting for 
a stochastic "annealed" growth probability on the crack perimeter. 



2.6.2 Crack roughness and fracture energy. The reasons why the rough- 
ening of cracks is interesting are multiple. First of all, as discussed below, 
there is the usual interest in possible universality classes depending on di- 
mension, fracture mode, and material (disorder, average response). Second, 
it has also an engineering implication. Consider a crack which "starts flat" 
and then undergoes a process which indicates a scaling behavior similar to 
Eq. (|25() . An infinitesimal crack opening will now cost more surface energy 
than in the usual case, and thus the Griffith argument (Eq. (|15() ) is modified 
in a way that accounts for the roughness exponent the maximal extent of 
fluctuations ^max, and possibly also the modified elastic energy due to the 
rough crack. In the case of anomalous scaling (see Eq. (|28() again), the extra 
scale also would enter the argument, further changing the amount of surface 
energy needed for infinitesimal crack advancement. This implies a non-trivial 
"R-curve", or crack resistance, since the energy consumed via forming new 
crack surface depends not only on the full set of exponents but again on the 
distance propagated [46,47]. 

An expression can be written down using solely geometric arguments and 
the exponents as given below: 

AM/,y)_A|^^^^^ if l»ay) 

where ^(y) = By^/^ is the distance to the initial notch y. There are two 
different regimes: for I ^ Hv)-, Ah{l,y) ~ y''^^, and for / <C C{y), which 
is characterized by the exponent (C — Cioc)/^- Important here is presence of 
anomalous roughening so that Ah{l,y ^ ysat) — Al'''"" L''~^'°'=. 
Since the effective area from which surface energy follows is 

G 6Ap = 27 6Ar (31) 

where 5Ar is the real area increment and 5Ap its projection on the fracture 
mean plane, the growth of roughness implies 



GniAa « Aasat) ^ 27^1 + (^^^^) Aa^C^-^'-^^ 
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and the resistance to fracture growth becomes 



For rough cracks, this is the simplest analytical attempt to consider size-effects 
(take l/L constant) and the effect of exponents describing self-affine roughness 
on fracture energy. 

Many authors have argued for the presence of "fractal" effects in the energy 
consumption (see e.g. [40]), while it seems clear that the relevant framework 
for calculating the true surface energy (7^) for the energy balance in the Grif- 
fith argument is a self-affine one [48, 49] (here we neglect the complications 
of multiscaling or anomalous scaling for simplicity). This can be argued for 
simply since the experimental situation seems clear-cut: usually cracks are self- 
affine (if not flat, but see also Ref. [50]). The attraction of (self-affine) fractal 
scaling ideas here is that they allow to understand the energy consumption 
mechanism(s) simply via an interpretation of the roughness exponent ^ or the 
fractal dimension Df. The models to be discussed in later sections are a tool- 
box for testing the influence of disorder on strength, since these models lead 
to rough, self-affine cracks. If one allows for mesoscopic plasticity, one would 
still expect the same kind of problems, but further theoretical questions are 
expected. 

3 Experimental background 

To give a reference point for the theoretical and numerical ideas discussed later, 
we next overview shortly some relevant experimental quantities, focusing on 
fracture properties for which disorder and stochasticity are crucial. When nec- 
essary we refer to the concepts presented in the Introduction, or point forward 
to the theoretical state-of-the-art presented in later sections. We would like to 
emphasize rather generally that the statistical aspects of fracture processes - 
scaling, exponents, distributions - are not understood even in the sense of a 
synthesis of experiments, with the main exception of the roughness behavior. 
This is by force reflected in the discussion that follows. 

The first of the issues is strength, both in terms of its probability distribu- 
tions and the scaling with sample size. As noted in the previous section, from 
the theory viewpoint, the statistics of strength depend crucially on the pres- 
ence of disorder. Next, we consider crack roughness, which has an application 
to the question of the origins of the fracture toughness. As pointed out already 
in Sec. 12.6.^ this is due to at least two reasons. The energy consumption that 
takes place when a crack increases in size is due to the product of surface 
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energy and the new area. Both can be affected by disorder, as in the case of 
Griffith-hke energy arguments and in the idea of crack surface area increasing 
with scale for self-affine cracks. In the case of fracture strength, disorder can be 
easily reinterpreted as a "flaw population", whereas for cracks that roughen, 
disorder plays the role of "noise" during crack dynamics. This noise can either 
be "quenched", i.e., exists in the sample from the beginning in many disguises, 
or be created dynamically via for instance void-formation. 

A very important topic from our point of view is represented by acoustic 
emission (AE) in fracture experiments, a nice example of crackling noise. The 
released elastic energy in fracture is a revealing diagnostic for the temporal 
aspects of fracture processes, and is used in applications as a monitoring tool 
for damage assessment. In addition, it has been long since realized that the AE 
signal appears to be critical: the energy of AE events often exhibits a power-law 
-like probability distribution. In contrast to thermal imaging or displacement 
field analysis, AE is not limited to sample surfaces and is sensitive to even 
very minute events. The two fundamental questions that concern us are: i) 
how to connect the AE behavior to the other fracture properties, and ii) how 
to understand theoretically the AE statistics. 

In the last part of this section, we briefly outline phenomena that exhibit 
time-dependent response, and often go beyond the scope of brittle fracture. 
This concerns the role of microscopic but collective effects on crack growth re- 
sistance ("R-curve") due to plasticity, and slow fracture processes with history 
effects such as creep and fatigue failure. We mostly omit discussing the spatial 
aspects of plastic deformation such as shear bands that bear some similarity 
to damage localization and fracture development in brittle and quasi-brittle 
fracture. 

3.1 Strength distributions and size effects 

Measuring the scale-dependent properties of strength is difficult mainly be- 
cause of sample-to-sample fluctuations. In addition, samples have free sur- 
faces, and though one often compensates by adjusting the test geometry for 
the resulting free boundary effects (Poisson contraction), materials often dis- 
play surface defects, giving rise to additional disorder. We can thus expect 
that the strength distribution will be a superposition of cases where surface- 
induced defects lead to failure and cases where bulk defects rule the process. 
Obviously, separating the two mechanisms can be difficult. In the classical ap- 
proach, the sample strength and the size effect is measured with a prepared 
notch. To obtain reliable results, one should have access to a large enough 
range of sample sizes. 

Fig. [7| presents experimental data from brittle ceramics, wherein the data 
are plotted according to different extremal statistics distributions [51]. In this 
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particular case, the Gumbel form fits the data better than the Weibuh distri- 
bution [51]. The latter is traditionally used in the engineering literature in the 
analysis of fracture statistics and size effects in strength [52,53]. One should 
pay attention to the fact that a Weibull scaling can be ascertained more re- 
liably by tests with varying sample sizes. One should note (see also Sec. 12.4(1 
that there is a priori no reason why a material should fulfill the prerequisites of 
belonging to any of the classes following from extremal statistics, be it Weibull 
or Gumbel. 

Various recent attempts to characterize the strength distribution in materi- 
als with intrinsic disorder such as paper [54,55], concrete [56] and others [57,58] 
illustrate clearly the inherent difficulties in interpreting the data. Fig.|Hldemon- 
strates this for concrete samples with a 1:32 size range. In this particular case, 
in addition to the inherent disorder and randomness, one has to deal with sys- 
tematic stress gradients present in these samples. This appears to be a generic 
problem that persists among various experimental studies (see e.g. [59]). 

In the engineering mechanics literature, the scaling of strength with sample 
size is usually discussed with the aid of three-point bending or tensile tests, 
done on samples in which the ratio ao of notch size to sample size is constant. 
Clearly, in the absence of any other effects, the square-root -like enhancement 
K ~ y/oQ for the stress intensity factor will eventually win the competition 
over the disorder, so that the strength will decay as expected by an inverse 
square-root law: cTc ~ l/y/oQ. For small enough oq the situation is, however, 
not so simple since the stress enhancement due to the notch competes with 
the disorder. In this respect. Fig. 1^1 demonstrates an example of data from 
concrete, where one can observe clear deviation from the ideal behavior. 

In the presence of disorder and/or plasticity the scaling with a constant 
aspect ratio may be complicated by one issue. Using engineering mechanics 
language, the fracture process zone (FPZ) size £,fpz may not be constant as a 
function of and L. Thus the energy balance (Griffith's) argument develops a 
size-dependence, which can be understood considering that the energy needed 
for crack growth changes with the volume of the FPZ. Such complications 
should depend on the details of the processes that make the FPZ change 
size. There is an extensive literature on the connection of microscopic and 
mesoscopic mechanisms to the size-effect, analyzing the possible behaviors of 
^FPZ for various materials. Many of these exhibit "strain-softening", but other 
processes as crack bridging by fibers or precipitates etc. are surprisingly general 
as well from composites to concrete [60-75]. This kind of analysis often also 
takes into account the role of free surfaces, which modify the FPZ geometry 
if the crack is close to a free surface. From our viewpoint, one should note 
that the presence of disorder can simply make small enough cracks irrelevant 
[76], masking the presence of such defects. Simply put, the randomness can 
correspond to an effective flaw size a^, so that for oq < a flaw is so small 
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as to not change the general behavior. 

In conclusion, the issue of fracture strength from the statistical physics view- 
point appears to still lack, to a large degree, physics-oriented experiments, 
which would help to connect scaling theories to reality. It would for instance 
be of immense interest to have empirical characterizations of the "disorder" 
(as a microcrack population) from sample to sample, and a concomitant mea- 
surement of the sample strengths. 




Figure 7. Examples of trial scaling plots for strength distributions of porous, brittle ceramics 
(from Ref. [51]). By plotting the distributions in various scales one can try to establish whether the 
shape resembles the Gumbel or WeibuU form the best. 
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Figure 8. The size-dependence of concrete for two different ambient conditions (adapted 
from [56]' ' ' . . „ ,, . ,. „ . „ ,^ ^.^^ 

D (for D ;ture. 




Figure 9. The typical size-effect from various experiments on concrete. It is clear that the data for 
fioxural strength, fr, exhibits a cross-over to a scaling as a function of D, that may (perhaps) 
follow the WeibuU size effect for large D values [74]. 

3.2 Rough cracks 

There are many different experimental signs of self-affine scaling for crack sur- 
faces in various scenarios reminiscent of those overviewed in Sec. 12.6.11 An 
excellent review by Elisabeth Bouchaud summarizes the major results up to 
1998 [77]. In general, one has two main issues: first, understanding the seal- 
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ing, depending on the details such as spatial dimension, material character- 
istics and loading state. In the language of self-affine geometry, this implies 
determining experimentally the scaling exponents, such as the /3, z, and C 
of the simplest Family- Vicsek scaling. Ultimately, a determination of the be- 
havior in terms of the exponents and scaling functions would point out the 
right stochastic equation that governs the process, including the noise that 
is present. Second, whatever the truth be about universality, one can search 
for the mechanisms that lead to the formation of rough interfaces, and try to 
understand the consequences of self-affine crack surfaces. In other words, re- 
gardless of whether one can write down an equation describing the dynamics of 
a crack front, the process is intimately coupled to other physical characteristics 
as acoustic emission, and damage mechanics. We classify various experimental 
results according to the dimensionality of the scenario described in Sec. I2.6.1L 
The conceptually simpler case of crack roughening is in general provided by 
case (i), where a the final crack path is given by a one-dimensional curve h{x) 
in two dimensions. Two typical examples are shown in Figures nUl and [TTl from 
paper fracture. These illustrate the diffusive nature of fracture, such that it is 
clear that the asymptotic scaling of h{x) can only be attained beyond some 
microscopic lengthscale in x and h. That is, this might imply that h{x) 
exhibits rms fluctuations (local width) with w{l) ~ only for I 3> ^m- An 
example is provided of such behavior in Fig. 1111 




Figure 10. The diffuse crack surface in paper (a two-dimensional material) and the thresholded 
"interface", shifted down for clarity (courtesy of L. Salminen, HUT). Here, the issue concerns how 
to define an unique h{x) in the case the damage field is continuous (one can see individual fibers in 

the image on the other hand). 

More quantitatively, the existing results point out towards a 2D (local) 
roughness exponent in the range C — 0.6 — 0.7 [78-80]. A number of ob- 
servations can be made: the measured value is higher than a random walk or 
a Brownian motion value (rw = 1/2, and ( > 1/2 implies interesting posi- 
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Figure 11. An example of a 2D crack-line from a paper web, of width 6500 mm [78]. The sample 
has an effective exponent of about 0.6 as shown by various measures, but only above a scale of 
about 2-3 mm (oc 30 pixels in the lateral direction). 



tive correlations for the motion of the crack tip. The value is also close to the 
numerical predictions of various models, both theoretical ones and the RFM 
as discussed in Sec. IH.4I However, even for ordinary, industrial paper itself 
there are numerous values available that are significantly higher than say 0.7 
(examples are found in Ref. [81]). It is not known at this time whether this 
reflects the difficulty of experimentally measuring or the fact that the result 
is not "universal" but depends on material parameters as for instance ductil- 
ity and anisotropy. In addition Bouchbinder et al. have indicated a scenario in 
which the h{x) has more complicated structure, exhibiting multifractal behav- 
ior (implying a gth-order dependent for the various q'-momenta of roughness 
measures, see Sec. 12. 6. T|) [82]. An alternative explanation is simply a crossover 
below ^rn as suggested above [78]. Note that the values of ( can not simply 
be interpreted as fractional Brownian motion, since history effects in h{x) are 
more complicated. 

Complicating our understanding of the problem is the fact that the role of 
kinetics (velocity, dependence on the loading scenario) has not been studied 
systematically. When experiments are started from a notch, the useful con- 
ceptualization is to consider the h{x) to be a trail of a particle, the crack 
tip. Then, there are two issues of immediate relevance: the role of the average 
crack tip velocity, and the translational invariance of the profile. The former 
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is important since various kinetic effects migiit be crucial, even for slow failure 
as in, for example, creep or fatigue fracture. The crack tip velocity is even 
more relevant for fast fracture where the timescales of sound waves are not 
well separated from the velocity. In analogy to the discussion in Sec. 12. 6. Ti the 
crack dynamics is expected to follow a scaling that depends on the distance 
to the crack tip, so that translational invariance is broken for x < ^, before 
eventual saturation. This caveat is in particular relevant if one considers a case 
where there is a FPZ whose size, and thus changes with crack length /. 

When a crack front propagates between two elastic blocks, the dynamics is 
enforced to be quasi-two-dimensional and one can again define an interface 
h{x) separating the intact from the failed region. The theoretical models for 
this experimental setup, referred as case (ii) in Sec. I2.H.H are dealt with sep- 
arately in Sec. The interest here lies with the roughness exhibited by h 
and the associated dynamics as h changes slowly with some average veloc- 
ity. Knut Maloy and collaborators have performed experiments on a couple of 
joined Plexiglas plates, previously sandblasted to induce a randomness. The 
dynamics of the crack is followed using a camera. As demonstrated in Fig. I12| 
the crack front moves by avalanches, whose areas are distributed according 
to a power law distribution P{A) ~ A~^'^ (Fig- The authors have also 
estimated the roughness exponent using different methods and the result is 
C ~ 0.6 [83-87]. Notice that this value is similar to the one measured in case 
(i) but the geometry and the elastic interactions are completely different. Thus 
there is no apparent reason to identify the two values. 

Most studies performed in the past, starting from the first pioneering work of 
Mandelbrot et al. [88], refer to three dimensional fracture, case (iii) according 
to Sec. 12.6.11 The evidence collected during the first decade of experiments on 
several materials suggested the presence of a universal roughness exponent in 
the range ^ — 0.75 — 0.85 [77]. The scaling regime over which this result is 
obtained is quite impressive, spanning five decades as it is shown in Fig. 1141 
In some cases, namely metals, it has been shown that a different exponent 
(i.e. C — 0-4 — 0.5) describes the roughness at small scales, with a crossover 
to the large-scale value. Fig. displays an example of both the short-scale 
and large-scale regimes in a typical (fatigue) experiment. This crossover has 
been initially explained as a dynamic effect: the small scale exponent would 
reflect the quasistatic limit and the large scale exponent the effect of finite 
velocities [77]. Notice, however, that the short-scale value is sometimes not 
observed, even at very low velocities as for silica glass [89], while in other 
cases, namely in granite and sandstones, the large-scale exponent is not seen 
even at high velocities [90]. 

Since the crack surface is separated from the intact part of a sample by a line 
that can fluctuate both in the out-of-plane and in-plane directions, one can 
in principle measure two exponents C|| and C± (see Sec. I2.6.T|) . This was first 
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Figure 12. Examples of crack fronts in a quasi-two-dimensional setup where a rough interface 
between two elastic plates is used as the substrate for the crack line. Notice the avalanche-like 
localized advancement of the line between subsequent frames. These can be used to measure locally 
the roughness exponent via the relation A.h ^ /'> where A.h is the average displacement, and I is 
the lateral extent of the region which moved [87] . 




Figure 13. The avalanche size distribution P{A) in the Plexiglas experiment of Ref. [87]. Notice 
the independence of the scaling exponent of the actual average crack velocity. 
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Figure 14. Roughness scaling measured in a metallic alloy broken in tension. Zmax{r), the 
maximum h-value difference inside a window of size r shows scaling over five decades with an 
exponent of ^ ~ 0.78. The inset shows the overall between two different measuring techniques (from 

Ref. [77]). 



attempted via elaborate experiments by stopping the experiment and cleaving 
the sample ( [77]). Very recent research has brought new light into these issues. 
Ponson, Bouchaud, and co-workers have tested the line propagation scenario 
on four materials with very different microstructural and brittleness/ductility 
properties [89,92]: silica glass, an aluminum alloy, mortar, and wood. Fig. [T^ 
shows an example of the correlation functions and the appropriate exponents, 
that seem to be valid for all these cases, and independent of the detailed mode 
of fracture (tensile, two-point bending, fatigue etc.). The parallel scaling is 
in agreement with earlier ^-results. It is noteworthy that in both cases one 
has the same dynamical exponent z, so that the two-dimensional correlation 
functions (to use the notation of Ponson et al., where x is the direction of 
propagation) would follow the scaling Ansatz. Ponson et al. find that z ~ 1.25 
which together with the Family- Vicsek relation then sets the three exponents 
(see Fig. HT)) 

Ah{Az,Ax) = Ax^f{Az/Ax^/^) 

, „, , Jl ifu«l (^^^ 

where /W~|^Cifu»i 

Thus, in this framework, the direction of crack propagation has the role of the 
time axis and is thus described by a growth exponent (3 = (±. This is because 
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Figure 15. Zmaxix) I ^f^-, the maximum fe-valuc difference inside a window of size r/^c for two 
materials. Here is a cross-over scale that can be determined separately and may be thought of as 
the FPZ scale. As can be seen for both cases and for various velocities of crack growth, there seems 

to be two different regimes of roughness with C, ~ 0.5 at small scales and C, ^ 0.8 at larger scales 

( [91]). 

the fracture surface results indeed from the trail left by the crack front and the 
dynamic scaling properties of the roughness are encoded in the perpendicular 
component. 

The role of microstructural disorder for fracture roughness is still debated. 
Recent experiments point out the difference of intergranular-transgranular 
fracture [93], and in many materials crack bridging should play a role in 
roughening. Some steps towards a more complete understanding come from 
the seemingly irrelevant role of crack growth velocity [94], and the fact that 
the same values for the roughness exponent can be measured from a very brit- 
tle material such as amorphous glass [95, 96] where the role of plasticity is 
negligible, and Cfpz being of the order of a hundred nanometers. Thus the 
roughness may originate from void formation, a process that has not been 
studied much from the viewpoint of statistical fracture. Bouchaud et al. [97] 
noted that in intergranular fracture, one sees similar geometric phenomena 
as in tree-like directed percolation clusters of statistical physics. The ideas of 
branching critical cracks have not been studied systematically since then. One 
should bear in mind that it is only the 2d projection of the 3d crack structure 
that is directly analogous to such a "percolation cluster". 
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Ax or Az (|im) 

Figure 16. ID height-height correlation functions for /ly and h^, in the notation of Sec. l2.6.'T1 for 
an aluminum alloy fracture surface [89]. The apparent scaling exponents are f = (Jy ~ 0.75 and 

/3 = C± ^ 0.58. 

We should also mention that there are signs of anomalous scaling from crack 
propagation experiments in three dimensions [33,98-100]. Thus the commonly 
measured local roughness at a distance I does not obey Family- Vicsek scaling, 
but we would need an additional independent exponent to describe the cou- 
pling to sample size, L, as discussed in Sec. 12.6.11 The origin of these scaling 
laws still remain to be understood, possibly in the framework of line depin- 
ning -models. Note that Bouchbinder et al. have recently suggested that the 
actual scaling picture is not simply self-affine, but one would need to resort to 
structure functions as in turbulence [101]. 

Finally, Fig. 1181 demonstrates how the crack growth resistance (as measured 
by the energy release rate Gptc) depends on sample size L in 3D experiments 
on two kinds of wood [46,47]. This can be compared to the theoretical scaling 
given in Sec. 12.6.11 This implies, that the self-affine picture of crack surface 
geometry would be directly useful in engineering applications if a mathematical 
relationship can be established between Gjic (or K) and C and the material 
in question. One should keep in mind two issues: the complete absence of a 
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Figure 17. 2D height-height correlation functions A.h/\^(h.x) for various Ax, Az for a fracture 
surface of an aluminum alloy ( [89]). The collapse uses the form suggested by Ponson et al., Eq. 1321 



theory to explain the value of C and the anomalous scaling picture, and the 
fact that quite possibly the microscopic mechanism that makes Grc increase 
with L also has a role to play in setting the actual value of C,. An engineer's 
way to explain this is that it is simply the prefactor that matters to first order, 
in K, not in particular if C, would prove to be universal. 

To conclude, we note that fractography is becoming a mature sub-branch of 
fracture mechanics, trying to establish links between the qualitative processes 
and the phenomenology of surface geometry on one hand, and on the other 
hand between crack topography and materials parameters such as fracture 
toughness and orientational effects (a good reference is the book by Hull [102]). 
There has been little effort to couple the statistical signatures of other quan- 
tities to those of crack surface roughening, however, regardless of whether one 
starts from the statistical physics or materials science viewpoints. 
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Figure 18. Size effect on critical energy release rate Guc for two kinds of samples: pine (a), and 
spruce (b). The expected slopes are using the associated roughness exponents C, — Cioc = 0.47 it 0.17 
for pine (a), and 0.73 ± 0.17 for spruce (b) (from [47]). 



3.3 Acoustic emission and avalanches 

The microscopic processes that take place in a fracturing system with one 
or with many cracks are difficult to follow experimentally. Some of the pos- 
sibilities available are thermal imaging, electromagnetic diagnostics, surface 
strain measurements including digital interferometry, elastic modulus diag- 
nostics [103], and three-dimensional tomography [104]. All of these have their 
disadvantages including lacking spatial and temporal resolution, inability to 
access bulk dynamics, and the presence of several noise sources that make 
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interpreting data difficult. 

Here we concentrate on perhaps tlie best way to characterize microfractur- 
ing dynamics, acoustic emission (AE). This is an example of "crackling noise" 
which has gained considerable attention in other contexts, from the view- 
point of statistical physics. Acoustic emission arises in brittle fracture due to 
the release of elastic energy, but can also originate from internal friction or 
dislocation motion, to name a few possible other causes. In the geophysics lit- 
erature [105] in particular it has been discussed in the context of scaling and 
patterning of fault systems and fracture networks [106-108], and the growth 
laws of the same [109,110]. In analogy to the simple statistical models discussed 
in this review the simplest theoretical interpretations concern the relation of 
released elastic energy to the detected AE. However, in most cases the models 
at hand do not take into account the discrete statistical and noisy nature of 
the AE data which is, here as well as in the context of earthquakes, really 
important. 

A typical AE experiment is based on detecting acoustic activity via piezo- 
electric sensors that convert the sound into an electrical signal. Various levels 
of sophistication exist, that allow to distinguish between exact waveforms and 
to triangulate the wavesource(s). After noise thresholding, the essential point 
is the existence of well-separated "events" such as depicted in Fig. It shows 
two subsequent events from a tensile test on a paper sample, and demonstrates 
the most important quantities that relate to such a signal: quiet times, event 
durations and event energies. In a slow fracture experiment one has a separa- 
tion of timescales, except perhaps close to the final catastrophic breakdown. 
The events as detected by an experimental system are difficult to interpret 
because of attenuation [111,112], measurement system response, and disper- 
sion including reflections from (internal) surfaces. Hence, only the simplest 
quantities such as energy, amplitude, and silent intervals are relatively easy 
to analyze. Nevertheless there is a considerable body of work on the relation 
of AE to microscopic fracture mechanisms and damage accumulation versus 
material composition, including for instance concrete and metal-matrix and 
fiber composites where the contrast between the fibers and the matrix, allow 
such a distinction [113-119]. 

Figure QUI depicts two tensile tests for paper samples, illustrating the typical 
dynamics of AE in a test. The interesting questions are, quite similar to crack 
roughness, how does the statistics depend on rate of loading, dimensionality, 
material, and fracture mode? At the simplest level the integrated acoustic 
energy should tell about the loss of elastic energy content, and thus be pro- 
portional to (1 — D) (in terms of a scalar damage variable) if one assumes 
that attenuation and acoustic coupling do not change as failure is approached. 
From a variety of experiments that include a variety of loading conditions in- 
cluding tension, compression, and shear, and in addition under creep loading 
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conditions, and different materials it has been found that the silent times r 
and the acoustic event energies E follow analogies of earthquake statistical 
laws. For the energy, one has 

P{E) ~ E-^, (33) 

and for the intervals 

P(r) ~ (34) 

where the exponents a and /? are akin to those of the Gutenberg- Richter and 
Omori laws of earthquakes. Typical values for a ~ 1 — 1.5 and /3 ~ 1 — 1.3 
are then a target for theoretical models [120-125]. An example of the energy 
statistics is given in Fig. 1211 for paper as an example of a two-dimensional test 
material. In an in-plane experiment, /? ~ 1.8 has been found [126] which is 
intriguingly close enough to the avalanche exponent of crack lines in Plexiglas 
(see Fig. ESI [87]). 

The example of Fig. QUI hinted already about one particular aspect of typical 
fracture AE: the lack of time-translational invariance since e.g. the activity 
gets stronger as the cTc is approached (though in special cases, as in "peeling" 
[126] this can be circumvented). This presents both a challenge in order to 
understand the statistics - recall that often one does not have access to large 
quantities of data - and implies that the AE timeseries contains clues about 
the microscopic dynamics. 

Fig. 122 (ref. [127]) demonstrates how the amplitude statistics changes as 
one follows a fracture experiment along the stress-strain curve. Here, the "b- 
value" or the analogous /3-exponent becomes smaller signaling that there are 
larger and larger energy releases as ac is approached. The b- value is one of 
the targets of theoretical analysis based on averaged damage and crack-growth 
laws, and is based on a "iCe//", which would predict energy release rates as 
the catastrophic fracture gets closer. These models can of course be tuned so 
as to match with the general trends of detecting larger amplitudes as the final 
fracture is approached. Still, the agreement is qualitative and the possible 
"universality" of the signatures remains unclear [128-134]. Similar b- value 
studies have outlined e.g. the role of pre-existing cracks in statistics [135], 
followed damage [136] and discussed the crack velocity's role [137]. 

The changes in behavior as sample failure is approached and the concomi- 
tant changes in data can be discussed in a variety of ways. One way to consider 
the matter is to try to do actual localization experiments. These often suffer 
from the fact that dead-times of the apparatus and uncertainties in localization 
make it difficult to follow all the developments (in contrast to plastic defor- 
mation [138], or steady-state fracture where the possibility of a quasi-steady 
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state exists [126]). Nevertheless, it can be convincingly be demonstrated how 
damage localizes to the final process zone, or to the proximity of an advancing 
crack in the final stages of stable crack growth. 

One can also check when or up to which point the assumption of slow failure 
is valid, in that the events remain uncorrelated on the fast timescale of inverse 
sound velocity [139]. The spatial clustering of events has been extensively dis- 
cussed in rock mechanics [140-143] and correlation analysis is often performed 
on the AE timeseries [144,145]. The latter presents a technical problem due 
to the non-stationarity of the time-series, as the AE event density often most 
varies during the experiment, and thus one should apply caution. In the spa- 
tial case, one can also do entropy analysis to distinguish between random, 
Poissonian statistics and real clustering [122]. This is demonstrated in Fig. 
1241 where one can see the development of correlations in the locations of AE 
events "in the plane", when a planar sample is subjected to an effective tensile 
load. The final failure surface is circular, as can be seen in Fig. 124b ,. The same 
samples (of wood) also allow to look at the clustering of the event locations 
via "entropy analysis" (Fig. Elb) in analogy to what we discuss about damage 
development in the numerical models below. An open question is, whether one 
can by spatial analysis say something about the dynamics inside and forma- 
tion of the FPZ [146] and determine the "Representative Volume Element" 
size (see e.g. [147]). 

Another consideration is to link the AE activity to final crack formation and 
properties, but unfortunately so far this has not been attempted (in paper, 
one can relate the jerky advancement of a crack tip - which is easy to demon- 
strate, see Ref. [148] - to AE activity, though). Instead, several authors have 
considered the AE energy release rate, or the integrated energy as a simplest 
way to see the approach to final failure. This naturally tells about the corre- 
lations in the damage development, and the nature of the final crack growth. 
Conflicting evidence exists, and Figure EHl shows that in some cases experi- 
mental data can be used to justify the possibility of a "critical divergence", in 
that the data seems to fit a scenario where dE/dt ~ {t^ — t)~°', where tc is a 
failure mode -dependent critical time (e.g. equivalent to a critical strain). In 
some cases, this seems not to work so well (an exponential increase is rather 
indicated) [149], moreover such data has also been used to justify the presence 
of collective dynamics in the process via "log-periodic oscillations", [150,151]. 
Such collective phenomena, of unknown exact origin, would be of importance 
since they would allow to consider AE as a diagnostic for sample life-time or 
failure predictability. 
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Figure 19. A typical example of the crackling noise of acoustic emission (courtesy of L. Salminen, 
HUT). The two events that are detected are often well-separated in time and are characterized by 
duration intervals t, and energy Ei. Note the presence of noise in the setup, which means one has 

to threshold the data to define events. 

3.4 Time- dependent fracture and plasticity 

Creep, fatigue and ductile failure extend the theoretical questions one meets 
within brittle or quasi-brittle fracture. Dislocations, activated dynamics, and 
frictional response (fiber pull-out in fiber composites, crack surface interac- 
tions) can be combined with the presence of "disorder" (the emphasis of this 
article) to repeat many of the relevant questions: how to describe the statis- 
tics of strength, how to understand the properties of crack surfaces, and how 
to interpret damage mechanics via acoustic emission. In fact in many of the 
experimental references cited above one can not so easily label the material 
or test at hand as purely brittle or "ductile". In some cases the physics is 
changed "relatively little" : the existence of ductility and plastic deformation 
is interpretable in terms of a FPZ size ^fpz inside of which the divergence of 
the stress-field of LEFM is cut off. 

The two classical problems of fatigue and creep should also be affected by 
disorder. Indeed, the "Paris' law" which states that the lifetime of a sample 
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Figure 20. Two examples of acoustic data series (crosses and circles) for two strain rates of paper 
samples in tensile testing. The stress-strain curves change with the strain rate (from Ref. [125]). 
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Figure 21. Three different energy distributions from a tensile AE test for paper (e = 1.0 [%/min]). 
The data is separated in addition to the total distribution for the post- and pre-maximum stress 

parts (from Ref. [125]). 
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Figure 22. "b-values" or the corresponding amplitude distributions for acoustic emission from rock 
samples. It is noteworthy how the exponent of the cumulative distribution becomes smaller (from 
"stage 1") as the maximum point along the stress-strain curve is approached [127]. 



in stress-cycled experiment, with a variation /S.K for the stress intensity factor 
scales as ~ {A.K)'^ where n is the Paris' law exponent, can be measured 
in various non-crystalline materials which are certainly inhomogeneous [152]. 
Likewise, the physics of creep in disordered materials as paper or composites 
is affected in some yet to be understood way by the randomness. This can 
be seen in the stick-slip -like advancement of cracks in creep [153] and in 
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the statistics of creep times and the related acoustic emission [154]. In these 
directions, much more careful experimental work is called for. 

In the case of pure dislocation plasticity, outside of "fracture" there are col- 
lective effects (see the review of Zaiser for an outline [155]), which often need 
to be treated with the tools of statistical mechanics. Examples of collective be- 
havior are the Portevin-LeChatelier effect [156], and the temporal and spatial 
fluctuations in dislocation dynamics [138,157]. One would like to understand 
the laws of damage accumulation and how to define a RVE or correlation 
length [134,158-162]. The plastic deformation develops in non-crystalline ma- 
terials patterns, whose statistical structure is largely not understood (see e.g.. 
Ref. [163]) though simple models have been devised [164]. These questions 
could also be studied again via statistical mechanics ideas and models. 




Figure 23. The cumulative energy E, normalized to Emax, as a function of the reduced control 
parameter (t — t)/T, plotted in the proximity of the failure in the case of a constant applied load 
(pressure). The circles are the averages, for 9 wood samples. The solid line is a trial fit 

E = Eq (-^-^) According to Ref. [124], 7 = 0.26, is independent of the "load". Note the range of 

E in which the scaling seems to apply. 



4 Statistical models of failure 

In the following we establish the rules of the game. We consider models with 
interacting elements, whose coarse-grained behavior captures the key features 
of fracture processes in the presence of randomness, that can be either annealed 
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Figure 24. a): the microfractures localize (in 2D) as the control parameter is increased. The final 
panel has all the located events from the five equidistant intervals in five subplots (from Ref. [122]). 
Note that any localization experiment typically misses some events that can not be placed with 
certainty, b): The entropy, which decays roughly linearly from one (corresponding to random 
locations) as the final fracture is approached for the pressure (control parameter) value of P/Pc = 1. 

(thermal) or quenched (frozen). Moreover one may consider models in which 
one incorporates a history effect. A typical case would be where the dam- 
age accumulated by individual elements would be a monotonically increasing 
stochastic function of its stress-history. The simple rules used are then in- 
tended to capture the behavior of materials on the mesoscopic level ranging 
from brittle to perfectly plastic to viscoelastic. Another extension concerns 
models with inertia, so that sound waves are included. We also briefly take up 
the issue of connecting these abstract models to the finite element simulations 
practiced in engineering fracture mechanics and to atomistic simulations. Note 
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Figure 25. Random thresholds fuse network with triangular lattice topology. Each of the bonds in 
the network is a fuse with unit electrical conductance and a breaking threshold randomly assigned 
based on a probability distribution (for example, a uniform distribution). The behavior of the fuse 
is linear up to the breaking threshold. Periodic boundary conditions are applied in the horizontal 
direction and a unit voltage difference, V = 1, is applied between the top and bottom of lattice 
system bus bars. As the current / flowing through the lattice network is increased, the fuses will 

burn out one by one. 

again that in this review we on purpose leave aside the well-estabhshed field 
of dielectric breakdown, in spite of its close connections to scalar fracture [22]. 



4.1 Random fuse networks: brittle and plastic 

The simplest truly interacting fracture models are variations of the scalar 
analogy of fracture captured by the random fuse model (RFM) [8, 12]. The 
equations governing the RFM can be considered as a discretization of the 
continuum Laplace equation 
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with appropriate boundary conditions (see Fig. l25() . This scalar electrical anal- 
ogy of elasticity is a simplification over the Lame equations (Eq. (fTU)) ) de- 
scribed in Sec. El and formally corresponds to an antiplanar shear deformation 
scenario. In the discrete version, the currents and voltages over the nodes Xij 
satisfy the Kirchhoff and Ohm's laws and depend on the local conductivities 
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cTjj for the links (or resistors between the nodes). The fuse is then simply an 
element with a property of irreversibility (see Fig. 126(1 : once a critical current 
(or similarly a voltage drop across the fuse, but if cjjj = 1 the two cases are 
perfectly equivalent) ic is reached, the conductivity goes to zero. Clearly, a 
model with such fuses is appropriate to describe the failure of very brittle me- 
dia. The failure becomes more gradual in the presence of disorder, since not 
all the links would fail at the same time. Although much of our focus is on 
the RFM, many of the issues that we consider are directly applicable to other 
discrete models as well. 

Quenched randomness is usually introduced in the RFM in different ways: 

(i) random thresholds, extracting the thresholds ic,ij of each link from a prob- 
ability distribution p{x); 

(ii) random dilution, removing a fraction p of the links at the beginning of the 
simulation; 

(iii) random conductivities, extracting the conductivities cjjj of each link from 
a probability distribution Pa{x). 

All these three methods allow to tune the amount of disorder present in the 
system by changing the distributions or the dilution ratio. 

The crossover from weak to strong disorder is mostly studied in the 
framework random thresholds (case 1) using either a uniform distribution in 
[1 — R,l + R] [165] or a power law distribution p{x) ~ x~^^^^^ in [0, 1] [166]. 
The first case interpolates between no disorder (R = 0) to strong disorder 
(R = 1), while the second case allows for extremely strong disorder, since 
for A — > oo the distribution becomes non-normalizable. The case R = 1 is 
probably the most studied in the literature. Alternatively it is also possible 
in principle to connect the disorder in the model to a realistic material mi- 
crostructure by prescribing the fuse properties based on the morphology of the 
material microstructure (see Sec. 14. 3|) . The effect of percolation type disorder 
(random dilution, case 2) on fracture will be explored in Sec. 15.41 The scenario 
of quenched random conductivities poses an interesting question on how dif- 
ferent disorder phases would need to be defined in the case of a RFM [167]. In 
general, analogous to elastic media with locally varying elastic constants [168] 
or simply percolative disorder, quenched conductivity variations create corre- 
lations among local currents i even in the linear IV-regime. The largest current 
or stress variations are related to specific, "funnel"-like defects (see e.g. [169]), 
which plays an important role in the weak disorder limit. Similarly, the use 
of irregular lattices is analogous to a percolation type disorder. These type 
of irregular lattice systems have been carried out using spring and beam-type 
lattice models [20,170]. 

The brittle RFM scenario is usually considered in the limit of time scale 
separation. This implies that the equilibration time of the currents is much 
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Figure 26. Examples of perfectly brittle and plastic fuses. A fuse undergoes an irreversible 
transformation from (V, I)c to {Vc, 0) in the case of a brittle fuse or keeps its current at the yield 
current [Ic or ly) in the case of a plastic fuse. In the brittle fuse case the response is linear up to 
{V,I)c- On the other hand, in the plastic fuse case, the unloading tangent modulus is same as the 
original tangent modulus, and in the fully unloaded state (zero current) there is a remnant "strain" 

or voltage in the bond. 

faster than the ramp- up rate of external potential or current. Moreover, any 
type of possible overshoots in the local currents it^ij are neglected (see Sec. 14.41 
for discussion on inertial effects). The dynamics is then defined by finding the 
quantity 

maxjj (36) 

at each timestep t. This is an example of extremal dynamics, similar to many 
toy-models that result in "avalanches" of activity. Here the physics ensues from 
the combination of local current enhancements and disorder, in the thresholds 
ic- The question is then, whether the latter plays any role or not, and, if it 
does, what are the consequences for a statistical and dynamical properties of 
fracture. Instead of successive bond failures of one at a time, for computational 
convenience, it is also possible to increase the external voltage in steps, and 
then let all the fuses for which the ratio of Eq. ()36() exceeds unity fail at once. 
This of course neglects the small adjustments of the static electric/elastic field 
that would take place in the usual one- by-one fuse failure dynamics. 

This discussion should at this point be made more precise by defining the 
pertinent quantities. Figure HZI shows a schematic of the I — V behavior of 
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a RFM. Note that the series of fuse failures defined by Eq. (|36jl is unique 
given the set of thresholds ic- The I — V characteristics depends, however, on 
the way the system is driven: two different I — V curves ensue depending on 
whether the system is driven by constantly increasing current or voltage. The 
figure implies that in the voltage driven scenario, there is the possibility of a 
weakening tail (the regime of / < Imax and V > V{Imax)), whereas in the 
current driven scenario, the curve stops at (V, I)max- Hence, depending on the 
loading scenario, there are two failure points. Under current control all the 
fuses in the post-peak regime (beyond the peak current Imax) break at once, 
whereas in the voltage case fuses in the post-peak regime fail gradually. The 
schematic presented in Fig. corresponds to a typical simulation in which 
failure of one bond at a time occurs. In order to obtain the average I — V 
response, a naive averaging of the I—V responses at constant number of broken 
bonds may not be the best option since such an averaging will smoothen out 
the maximum, exactly as if many such systems were in parallel. One way 
to overcome this problem would be to center first the voltage axis on the 
breakdown point V{Imax) of each realization and then average. Alternatively, 
an average of the envelopes of the I — V curves can be performed at constant 
voltage values (see Sec. 16. If) . 

An interesting variant of the RFM is given by the perfectly plastic one. The 
corresponding fuse response is pictured in Fig. 1261 a blown fuse, instead of 
failing, "yields" so that its current is set to a constant. In this case, the yield 
iteration proceeds until there is a connected yield path crossing the system, 
which then blocks any further increase. These systems can be simulated via 
a "tangent-algorithm" envisioned by Roux and Hansen. The crucial point is 
that after a fuse yields any further progress along the I — V curve can be 
considered in a state in which the yielded fuses are simply removed. Then, the 
I—V behavior will be composed of piece- wise linear parts. The important thing 
here is the it,ij/iy,ij ratio, and during each step one subtracts from iy.j the 
already existing current. Meanwhile permanent strain (or actually voltage) Vpi 
develops. The tangent modulus or conductivity stays still at the initial value, 
which means that if the system is allowed to unload to zero current then there 
will be a remnant voltage, Vpi^ in each of the plastic fuses. Consequently, local 
internal stresses develop in each of the fuse of the lattice system. All this 
rich behavior arises in a model with simple scalar plasticity and without any 
time-dependent deformation. 

Since the final yield current is connected to the presence of a yield surface, 
it is possible to understand its properties without resorting to actual RFM 
simulations. The path is given by the solution of the following optimization 
problem: minimize the value of the path-cost, over all possible configurations. 
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Figure 27. The schematic of a current-voltage curve using RFM. Two alternatives, current- and 
voltage-driven, are outlined. The schematic corresponds to a failure process in which one bond at a 

time fails. 

when the path energy is given by 

{<jy)L = ^ CTijepath- (37) 

path 

This is exactly the energy of the random bond Ising Hamiltonian at zero tem- 
perature (groundstate), so a number of properties are apparent (see [171] for 
a review). The yield surfaces in 2D and 3-D are self-affine, with the scaling ex- 
ponents C2D = 2/3 and CsD ~ 0.41. These can be affirmed by using a mapping 
(in 3-D) to the max-flow/min-cut problem of combinatorial optimization, and 
in 2-D the task maps to the directed polymer in a random medium (DPRM) 
partition function, which can be connected to the Kardar-Parisi-Zhang equa- 
tion (KPZ) of kinetic interface roughening [172]. The KPZ problem is exactly 
solvable, and thus the (20 follows. Further consequences for yield surfaces are 
the fact that the yield strength fluctuates, with the standard deviation scaling 
as L^, where 29 = 2^ + d — 3, and d is the dimension {d = 2 for 2D and 
d = 3 for 3D). Also, the exact nature of the probability distribution func- 
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tion (pdf) of fuse yield thresholds is irrelevant, as long as the thresholds are 
not too strongly spatially correlated, and the pdf does not have too fat a tail 
(power-law distribution, scaling as P{iy) ~ iy"" with a < 4). 

Finally, it should be noted that simulations of brittle-plastic networks of the 
kind depicted in Fig.ESlhave not been carried out so far. Such networks would 
present interesting problems to analyze from the theory viewpoint since the 
scaling properties of fracture surfaces, say C,2D-, are expected to be different. 
From a computational point of view, one problem is that unloading due to 
fuse failure causes yielded elements to return to elastic/linear behavior. Fur- 
thermore, returning to zero external load nevertheless creates internal stresses 
along with irreversible deformation [173]. 

It is interesting to mention here that there is a way to interpolate between 
fully brittle and plastic response. Instead of removing a fuse, it is possible to 
just reduce its conductivity (akin to softening) by Uij — > (1 — D)aij, where 
D refers to damage. This leads naturally to crack arrest [174], and thus to 
distributed damage. If one iterates the procedure ( [175]) so that continuous 
damage (Eq. ((201)) is created by several fracture events per bond, one can 
observe an effectively plastic I — V curve. Note that the tangent modulus is 
reduced, but there is no irreversible deformation. Models of this type have been 
found effective to model biological materials [176], but could also be useful for 
modeling softening of concrete. 

Finally we notice that there are some experimental studies of systems cor- 
responding to realizations of the RFM. Otomar et al. build a network of fuses 
similar to the numerical models considered in Sec. 101 [177], while a thermal fuse 
model has been worked out experimentally in Ref. [178]. Such attempts, while 
being laudatory, also highlight the difficulty of reaching substantial system 
sizes. Otomar et al. have studied the I — V curves and the dependence of the 
strength on the effective disorder, which was varied by putting together sam- 
ples differently. Their experimental results are all comparable to the numerical 
results presented in Sec. |01 



4.2 Tensorial models 

Following the general strategy of the RFM, several models with more realistic 
elastic interactions (see Sec. l2.1() have been proposed in the past. The simplest 
possibility is provided by central force systems, where nodes are connected by 
elastic springs [179-182]. The random spring model (RSM) is defined by the 
Hamiltonian 




(38) 
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where Ui is the displacement of node i and K the spring constant. The elastic 
equilibrium is obtained by minimizing Eq. H38|). and disorder is introduced in 
the standard way by either random dilution, random threshold, or random 
elasticity. Despite the fact that the RSM suffers from the caveats associated 
with the presence of soft modes, it is still considered quite often in the context 
of fracture. In terms of continuum elasticity, the discrete triangular central- 
force lattice model represents an isotropic elastic medium with a fixed Poisson's 
ratio of 1/3 in two-dimensions, and 1/4 in three-dimensions. 

A similar slightly more complicated model is the Born model, defined by the 
Hamiltonian 



H 



^ ^ \a{ui - tiJOy + f^{ui - Uj)\ . (39) 



2 



where (i, j) represents the nearest neighbors of the node i, and the summation 
is carried out over all the nodes in the lattice system. The quantities {ui — 
Uj)\\ and {ui — Uj)i. represent the relative displacement of the node j in the 
directions parallel and perpendicular to the bond (i,j) respectively. In this 
case, there is a primitive competition between stretching and bending the 
"bond" between two nodes on a lattice. Note that in the limit a = [5 the Born 
model reduces to the random spring network (if u is scalar, then it reduces to 
random fuse network), but otherwise it still lacks one fundamental property of 
real elastic systems that local rigid rotations cost energy (see Refs. [183,184]). 

The most used (and most physically correct) models in this respect are the 
so-called beam and bond-bending models. One may view this kind of a system 
as a lattice of massless particles interacting (or connected) by a beam to their 
nearest neighbors [185] . In the beam models, each beam is capable of sustaining 
longitudinal {F) and shear forces {S), and a bending moment (M). Following 
the notation of Ref. [186], one can define three effective parameters a = l/EA, 
b = l/GA, and c = /EI. For a more generic treatment of 2D and 3D beam 
models, refer to Ref. [187]. Using standard beam theory, it is possible to relate 
the beam displacements and rotations {x and y displacements and rotation) 
to the three forces F, S, and M using the parameters a, 6, and c. The out- 
of-plane deformations or "buckling" of 2D beam lattice models have however 
not been dealt thoroughly in the literature so far. In the literature, the beam 
models have often been used to describe the fracture of random fiber networks, 
as lattice models of composites and paper. In terms of continuum elasticity, 
the beam models correspond to discretization of Cosserat continuum elasticity. 
The bond-bending model, on the other hand, attempts to be a generalization 
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of the central force model [181,188,189], whose Hamiltonian is given by 



where K and B are extensional and rotational rigidities. In the above Hamil- 
tonian, the first term corresponds to the usual extensional energy between 
neighbors i and j, and the second term corresponds to the changes in the an- 
gle between two "neighboring bonds" emanating from i to j and j to k. This 
has of course the effect of making such angles rigid and thus resisting local 
rotations of bonds. The corresponding reactions of the bond bending model 
are longitudinal (F) and shear forces (S) that can be sustained by each of the 
elements of the lattice, and a twisting moment or torque (M) at each of the 
nodes of the lattice system. 

The inclusion of bond-bending terms into a triangular lattice system result 
in a conventional 2D isotropic elastic medium with varying Poisson's ratio 
values [190]. In particular, even the square lattice system with bond bending 
terms and next to nearest neighbor interactions is equivalent to a 2D isotropic 
elastic continuum [190]. On the other hand, the triangular beam model lattice 
system represents a Cosserat continuum. 

The failure condition of these bond-bending and beam lattice models de- 
pends on the combination of local torque, longitudinal and shear forces in 
beam models, or on central and angular forces in the bond bending model. In 
the beam case a typical fracture criterion is of the type 



where tp and are force and moment thresholds respectively, which can 
be chosen separately for each bond. Similar to the RFM, one can then tune 
the system by scaling the external load so that only one bond fails at each 
time. In addition, it is possible to control the relative roles of the two fracture 
mechanisms by either adjusting the magnitudes oi tp and tM or adjusting the 
beam elasticity (slenderness or relative extensional and bending rigidities of 
the beams). It should be noted that Eq. (|1T|1 resembles the von Mises yield 
criterion. In this respect it is worth pointing out once again that there are no 
studies of plasticity or elastic-plastic fracture using these vectorial formulations 
for interactions. 




(40) 




(41) 
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4.3 Discrete Lattice versus Finite Element Modeling of Fracture 

Traditionally both discrete lattice models and standard continuum finite el- 
ement models are often used for modeling progressive damage evolution in 
quasi-brittle materials, leading to macroscopic fracture or crack propagation. 
In the discrete lattice models, the continuum is approximated a priori by a 
system of discrete elements such as central-force spring and beam elements. 

With regard to fracture, the standard continuum based finite elements are 
especially suitable for modeling progressive damage evolution or fracture of 
homogeneous (at least at the length scale of interest) materials. Simulation 
of fracture in heterogeneous materials is however complicated by the presence 
of disorder, whose presence naturally leads to statistical distributions of fail- 
ure stresses, accumulated damage, acoustic activity, crack shapes and so on. 
Standard continuum constitutive equations are based on mean-field or homog- 
enization considerations, and hence do not capture the effect of fluctuations 
as has been briefly discussed in Sec. |21 A few of the recent continuum based 
investigations however have made progress in explicitly modeling the material 
disorder to account for its influence on dynamic crack propagation [191-194]. 

In the discrete lattice modeling of fracture, the disorder can be explicitly 
modeled. As discussed in Sections 14.11 and 14.21 the medium is described by 
a discrete set of elastic bonds with randomly distributed failure thresholds, 
or in the electrical analogy via the fuse model. In the elastic case one can 
resort according to the continuum elastic description desired to bond-bending 
models, the RSM, or the beam-type models. These have been used extensively 
in the mesoscopic modeling of progressive damage evolution in quasi-brittle 
materials. Sections 14.11 and 14.21 describe how statistical effects of fracture and 
disorder may be modeled using these discrete lattice systems. Applications of 
central-force spring elements to model fracture of homophase and heterophase 
composite materials is presented in Refs. [195-197]. Similarly, Refs. [16-19] 
present the application of beam lattice models for modeling fracture of concrete 
structures. 

On the other hand, in the continuum mechanics based finite element mod- 
els, the disorder is quite often implicitly modeled using the homogenization 
concepts over RVE and by degrading the overall constitutive response (stiff- 
ness) of the material based on effective stress and damage variable defini- 
tions [17,24,25]. In the continuum sense, the fracture of quasi-brittle materi- 
als is studied from two different perspectives, namely, the continuum damage 
mechanics and the continuum fracture mechanics. The continuum damage me- 
chanics deals with the study of crack formation and growth from an initially 
flaw or defect free structure. In particular, it describes the progressive loss 
of material integrity due to the propagation and coalescence of microcracks, 
microvoids, and similar defects. On the other hand, the continuum fracture 
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mechanics studies how a preexisting crack-hke flaw or defect would grow. 

Two continuum damage mechanics based approaches are commonly em- 
ployed at the macroscopic level for modeling damage evolution. In these ap- 
proaches, the bulk behavior of the material is described by either a phenomeno- 
logical model [17, 24, 25] or a homogenized micro-mechanics based model. As 
described in Sec. [2 damage is described by a scalar damage state variable 
D (or a generic tensorial D). The phenomenological models are based on the 
notion of an effective stress and the constitutive response of the bulk mate- 
rial is described by a phenomenological strain energy function. An evolution 
law for damage is prescribed along the gradient of a damage loading surface. 
Broadly speaking, the continuum damage mechanics formulation parallels that 
of standard plasticity theories except that the effects of strain softening must 
be considered in modeling damage evolution leading to fracture. 

Homogenized micro-mechanics based damage models [198-201] consider an 
averaged response of the material microstructural features. In these models, 
the macroscopic damage description is obtained by homogenizing the mate- 
rial response over a representative volume element (RVE) unit cell. The ho- 
mogenization of the material response, however, is valid only as long as the 
material is statistically homogeneous. While accounting for short-range local 
stress fluctuations and multiple crack interactions, the homogenized micro- 
mechanics models suffer from the same shortcomings as the phenomenological 
damage models. In particular, these formulations do not include long-range 
spatial correlations of cracks and stress fluctuations [198]. In the crack growth 
and coalescence stages of damage evolution, statistical fluctuations become 
dominant and homogenization based approaches are not applicable. 

It is useful again to compare the RVE size, ajivE, to the correlation length. 

In practice, one needs ^ <C ajivE for the homogenization to be accurate. 
Beyond the correlation length, ^, the material can be regarded as statistically 
local and homogeneous. With the increasing value of correlation length, i.e., 
as ^^^^ — > 1, finite size corrections are necessary for homogenization, and the 
material response changes from a local to nonlocal behavior [198]. Within the 
scale ^, statistical fluctuations are dominant, and hence homogenization is not 
valid. Close to macroscopic fracture, the extreme moments of statistical distri- 
bution of microcracks and the corresponding correlation lengths dominate the 
material response. Under these circumstances, the densities become irrelevant 
and the concept of RVE ceases to make sense. This drawback in homogenized 
methods has enabled the more recent investigations to explicitly model the 
microstructural disorder in finite element simulations [192-194]. 

An example of continuum damage mechanics based finite element model 
applied to the fracture of disordered materials such as concrete is the smeared 
crack approach [17], wherein the cracks are not explicitly modeled. Instead, 
the constitutive properties of the finite elements where a crack is supposed to 
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develop are degraded according to some continuum law. That is, instead of 
considering the crack explicitly as a discrete displacement jump, the effect of 
crack is smeared over all the Gauss points of the finite elements that contain the 
crack by degrading the corresponding material properties at the Gauss points. 
Many such finite element models based on continuum damage mechanics are in 
common use for modeling damage evolution including strain softening behavior 
of materials. 

On the other hand, in continuum fracture mechanics based finite clement 
models such as the the discrete crack method [17], an explicit modeling of the 
crack is considered using appropriate interface elements and crack-tip finite 
elements. The zero thickness interface elements are in general embedded be- 
tween the edges of two adjacent finite elements and/or along probable crack 
paths that are decided upon a priori. The constitutive behavior of an inter- 
face element is governed by a traction-displacement cohesive law, which is used 
to evaluate interface traction in terms of normal and tangential displacement 
jumps, or a non-dimensional displacement jump parameter which is zero at the 
traction free state and equals unity at interface failure. The progressive failure 
of these interface elements simulates crack propagation. When the remeshing is 
not performed during crack propagation, the crack path is restricted along the 
edges of the finite elements. This results in a finite element mesh-orientation 
dependency of the numerical solution. On the other hand, remeshing in each 
load step may be considered to alleviate the problem of crack path propa- 
gation on mesh layout. However, remeshing techniques [202], which model a 
discontinuity or crack by modifying the finite element mesh topology to ex- 
plicitly capture a discontinuity, are computationally expensive and not readily 
applicable for non-linear problems. Alternative fictitious crack model [203] ap- 
proaches based on ideas using Dugdale-Barenblatt plastic crack-tip zones are 
also in common use; however these approaches suffer from similar drawbacks. 

During strain localization and softening, the deformation concentrates in 
regions of finite size. The modeling of this localization process using continuum 
finite elements is fraught with serious difficulties, both from mathematical 
and numerical points of view. Mathematically, the onset of strain softening 
leads to ill-posedness of the rate boundary vahic problem of the continuum 
equations [204]. In particular, the governing set of partial differential equations 
in the standard rate-independent continuum lose their ellipticity for quasi- 
static loading conditions, while hyperbolicity may be lost locally for dynamic 
conditions [204]. In analytical solutions, strain softening is characterized by 
the appearance of infinite strains of a set of measure zero. When the problem 
is modeled and discretized with finite elements, the localization concentrates 
in the volume of material capable of representing the set of measure zero, 
which is a one element in one-dimensional problems or a one element wide line 
in two-dimensional problems. Prom a numerical perspective, this ill-posedness 
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manifests itself in severe mesh dependence of the finite element solution in 
terms of mesh size and alignment. Hence, the problem of mesh dependency is 
not intrinsically numerical but stems from the change of the character of the 
underlying differential equations of the continuum description. 

In order to allow for modeling of localization and strain softening in the 
continuum description, approaches based on strain localization limiters are 
proposed. The essential idea is to change the character of the equations so 
that the region of localization docs not generate to the set of measure zero 
but to the physical region that can be deduced from the experimental proce- 
dures. In the literature, several regularization techniques (localization limiters) 
based on spatial nonlocal continuum concepts [205] , higher-order gradient the- 
ories [206], micro-polar or Cosserat continuum descriptions [204], viscoplastic 
regularization methodologies [207] , and the regularized discontinuous displace- 
ment approximations [208] have been presented for capturing the localization 
phenomena. These approaches ensure that the numerical solutions remain 
meaningful in the sense that a finite width of the localization zone and fi- 
nite energy dissipation are obtained even when the localization zone is fully 
developed. The above localization limiters (regularization techniques) can be 
classified as integral, differential and rate limiters. Integral or non-local lim- 
iters are based on integral strain measures over a finite domain. Differential 
limiters include strain or stress measures with derivatives of order higher than 
one. In rate limiters, the time dependence of the strain softening is built into 
the constitutive equations. Although, each of the above methodologies has its 
own advantages and disadvantages, from an implementation point of view, 
differential and rate limiters are the realistic candidates for implementation 
into existing codes. Non-local limiters are more complex both theoretically 
and implementation-wise. 

In order to eliminate the effects of mesh size and alignment dependencies 
in the numerical simulation of fracture, finite element methods based on so- 
called embedded discontinuities [208, 209] within the finite elements are also 
used. These finite element formulations are based on enhanced assumed strain 
(EAS) formulations [210], and the displacement jump across a crack is in- 
corporated into a finite element as an incompatible mode in the strain field. 
Recently, extended finite element methods (XFEM) [211] based on the parti- 
tion of unity approach [212] are also employed for simulating the fracture of 
quasi-brittle materials. Although the XFEM has similarities in the modeling 
of cracks using discontinuous displacement interpolations within the elements, 
the XFEM is far superior to the discrete discontinuous displacements based 
finite element approach. To start with, XFEM is capable of representing ar- 
bitrary propagation of cracks within the elements with or without bisecting 
the elements completely [211]. For example, in the mesoscopic fracture simula- 
tions, this feature allows for the simulation of transgranular crack propagation 
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in addition to the usual intergranular crack propagation [213]. Furthermore, 
XFEM allows for the presence of multiple cracks with secondary branches in- 
cluding their mutual interaction thereby determining the change in the crack 
orientation. It is possible to enhance these discontinuous interpolations with a 
suitably chosen near crack-tip fields to achieve a superior coarse mesh accuracy 
in the simulations. Representing the cracks using a signed distance function 
based on level set methodology facilitates the evolution of cracks and hence 
alleviates the necessity to use automatic remeshing techniques and explicit 
representation of cracks [211]. 

Even with these vastly improved finite element methodologies and the re- 
cent trend in explicitly modeling the material disorder [191-194], fracture 
simulation of disordered quasi-brittle materials with many interacting microc- 
racks, and long-range spatial crack correlations and stress fluctuations poses a 
daunting task for current finite element methodologies. In particular, when the 
macroscopic fracture is preceded by many competing and interacting microc- 
racks, both mathematical and numerical issues arise in assessing which crack 
should propagate and by how much should it propagate [214]. In this sense, the 
rigor and accuracy of these homogenized finite element models based on effec- 
tive medium (or mean-field) continuum theories is limited to moderate damage 
levels, whereas those based on continuum fracture mechanics is limited to the 
regime where the behavior is dominated by a single (or few) dominating crack. 
Currently, the most fundamental challenge is to bridge the gap between the 
discrete and continuum descriptions of brittle fracture. Recent finite element 
studies [191-194] that explicitly include the material disorder effects are very 
promising in this direction. 



4.4 Dynamic effects 

4.4.1 Annealed disorder and other thermal effects. In reality, on scales 
that are supposed to be described by the quasi-static models introduced above, 
fracture processes may also be influenced by time-dependent details. A possi- 
bility that has been widely explored in the past is to move from a deterministic 
dynamics in a disordered system to a stochastic dynamics in an homogeneous 
system [215,216]. 

Time dependent randomness has been implemented in two different ways 
in the REM: the failure probability can be chosen to be a function of the 
instantaneous fuse current, or a functional (such as an integral) that depends 
on the current history. The first case corresponds to annealed disorder, where 
probability of failure of fuses is given by [21] 



Pij ~ i.H,ijY'^ 



(42) 
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properly normalized over all ij, such that r] measures the relative role of cur- 
rent enhancements. In the limit rj = 0, the model is equivalent to screened 
percolation. In the opposite case of r/ oo, the current enhancements al- 
ways dominate over any randomness (or disorder). It is commonly believed 
that quenched and annealed disorder lead to similar behavior in such a sce- 
nario. The simplest explanation is that one can always find a set of quenched 
thresholds such that it will exactly reproduce the same sequence of failures as 
any stochastically induced history [166]. The behavior of the RFM and RSM 
for various rj values has been explored both using simulations and mean-field 
theory [217-225]. A critical value seems to be r/c = 2 that arises due to a com- 
petition of damage accumulation and crack growth. Many issues still remain 
open, of those to be discussed in Sec. 101 the crack roughness is an example. 

In the second case (history dependent failure probability scenario), one may 
consider the case where the local current it^ij is followed as a function of time. 
For example, consider the integrated current, or some power thereof to a fixed 
threshold, and then use a rule similar to Eq. (|42|) . Under these circumstances, 
the scenario mimics a creep-like failure due to the external timescale that is 
now introduced by the integrated history. Consequently, it is most natural to 
study such models at constant current or voltage since a ramp-up in loading 
will only complicate the behavior further. 

Alternatively, the dynamics may be studied by considering the evolution of 
an extra state variable, T, defined in a manner similar to temperature [226]. 
The dynamics of T is supposed to obey a state equation given by 

CdtT = Ri^ - aT. (43) 

Now, if we suppose that a fuse fails whenever T > Tc, where the critical value 
Tc may be distrbuted uniformly, one again observes two regimes: percolation 
when 6^0 and strict brittle fuse networks when 6 — > oo [226]. For b ^ oo, the 
studied examples involve quenched disorder via local conductance variations. 
In the intermediate range of 6's, the networks develop interesting damage 
dynamics. 

4.4.2 Sound waves and viscoelasticity. The dynamics just discussed do 
not explicitly involve the two relevant time scales of loading and relaxation. 
The relaxation time scale in general encompasses two aspects; namely, the 
time scales that arise due to transient dynamics, i.e, as the microscopic dam- 
age takes place ("bonds break") the stress field deviates from its static or 
asymptotic stationary value, and the relaxation time scales inherently asso- 
ciated with the local dissipative mechanisms such as viscoelasticity. In the 
lattice models, these local dissipation mechanisms may be introduced locally 
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on an individual element level. 

The transient dynamic effects consider coarse-grained sound waves together 
with inertial effects. The perturbations induced by the dynamics into a quasi- 
static model can be considered in a time-dependent model through a dissipa- 
tion mechanism. As an example, let us consider a simple scalar 2D model 



This could be considered as a dynamical version of a RFM with the masses 
M and the damping F included, or a bare-bones simplification of a central- 
force one. Here the sum of site runs over the nearest neighbors (l,m). A 
constant loading rate (velocities) of V and —V is applied at the two opposite 
boundaries so that the boundaries move in opposite directions. 

The transient dynamics of the system are followed by numerically integrat- 
ing the Eq. ()44() . which is a damped wave equation. Each time a bond is 
stretched beyond its threshold Uc^ij, the bond is removed from the lattice sys- 
tem. In this scalar model one has only transverse waves, with the sound speed 
c = sjKj p. With such dynamics, one may monitor the local signal fFig. I4.4T^ 
from an "event" which does have some qualitative resemblance with real acous- 
tic events from experiments (Sec.|21). One of the complications associated with 
transient dynamics simulations is that wave interference and reflections signif- 
icantly disrupt the accuracy of numerical solution. Consequently, the damping 
coefficient F must be chosen such that the waves are damped strongly enough 
to avoid interference and reflections. This becomes particularly important in 
problems involving dynamic fragmentation. Another problem commonly asso- 
ciated with transient dynamics is also that close to surface boundaries incom- 
ing and outgoing reflections may interfere disrupting the solution. 

In the context of elastic waves excited by energy release, there is a multitude 
of possible mechanisms that are of interest: 

• local, transient stresses - due to a wave passing over a point x - may induce 
local failures, similar to a "beam" failure in a lattice beam model. 

• interference from several waves and from the descendant fronts of the same 
wave as it gets reflected by already existing microcracks [228]. The main 
point is that how these wave interferences amplify the dynamical, fluctuating 
stress (t(x, t) > o static beyond its static value. 

The consequences of dynamics in controlled (if not completely quasi-static) 
fracture would most likely be felt in the advancement of crack lines in 'id 
systems; analysis of the related effects indicates that this may induce crack 
front roughening beyond the otherwise at most logarithmic scaling [229,230]. 
In models with distributed damage, the consequences thereof would be rather 
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Figure 28. An example of a local amplitude at a site i. Clearly, individual events can be 
distinguished with an envelope that decays as a function of time. The lowest panel presents the 
integrated wave energy as a function of time. With linearly increasing applied strain, the increase 
in local amplitude as the failure point is approached is also visible in the time-series of "events" 
presented in the topmost panels, (see Ref. [227]) 



unclear. One must also recall the attenuation of waves from single source 
event with distance. It seems to us that it would be of interest to investigate 
systematically the role of deviations from the quasi-static limit, in particular 
as the failure point is approached. 

The lattice model discussed above is one of many examples used to study 
dynamical processes. Some of these we have encountered above in the context 
of DEM, the Discrete Element Method. There is a rather large literature on 
applying the DEM (also under the pseudonym Particle Method) to various 
materials problems, in particular to that of fracture behavior of concrete. 

A more physics-related application is the study of fast, dynamic crack propa- 
gation in 2d (for an overview, see [231]). Originally, the question arose how such 
simple models [232] would relate to the crack velocity vs. time, its asymptotic 
value in relation to the Rayleigh velocity, and the possible pattern formation. 
The models are studied either by taking the continuum limit, by simulations. 



61 



or by considering in detail the effects induced by the discreteness. The phe- 
nomenon of lattice trapping is an important consequence of such studies. 

Various authors have elaborated on the crack branching, in particular on 
the role of crack velocity oscillations, damping, and the exact microscopic 
force relation and so forth. There is however little overlap with statistical 
mechanics aspects such as the disorder in the samples, thermal effects and 
the thermodynamic aspects (where does the released energy go?). The same 
holds for various phase-field models of fracture [233 236]. In these models, 
one combines two fields: one for the elastic deformations such that a minimal 
choice for elasticity (e.g. in 2d) holds, and another one that allows for a order 
parameter p{x) that continuously interpolates between undamaged (p = 1) 
and failed {p = 0) states. The models are designed in order to simulate fast 
fracture propagation and the concomitant pattern formation. In the simplest 
form, phase-field models do have the problem that a crack is formed due to 
the presence of "sharp interfaces", with an intrinsic thickness or length-scale. 
These interfaces separate regions of the phase field p = 1 where the material 
is intact from completely broken ones (p = 0). The interface width is a fixed 
parameter which gives a minimum FPZ scale. 

There are no attempts to apply such models to situations with randomness, 
and perhaps one may imagine that there are difficulties in justifying physically 
the choices applied for the noise that one would add. For example, it is con- 
ceivable that one might start from a fourth-order Ginsburg-Landau -style free 
energy describing the intact and failed phases, and couple that to a field term 
that establishes a local preference ("toughness" or local strength). However, 
it is not clear how to adjust the noise strength in relation to the intrinsic FPZ 
size. 

One further possibility is to induce dissipation in lattice models such that the 
time relaxation effects are directly incorporated into the equations of motion. 
This can be done simply on the model level by applying rules from dashpot- 
style models of viscoelasticity [237-239]. One possibility is the Maxwell model 
in which the force fij acting through bond ij is given by 



dfij ^ _ 1 

dt dt --"^ 



- -fij (45) 



Here r is a phenomenological damping constant, and has the effect of relax- 
ing the local forces, fn is the contribution of a Hamiltonian, induced by the 
displacements of nodes. One could use any particular choice for it, scalar or 
vector. The system behavior is controlled by the ratio of two timescales: r 
and the strain-rate (or rate of loading). In a model that includes inertia ef- 
fects, there exists another time scale related to the dynamics of the system, 
and hence many kinds of behavior could be foreseen. Other possible examples 
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instead of Eq. (|45jl would be Kelvin and Voigt -types, which include damping 
in parallel to a spring. The corresponding response functions are of the form 

e = a/E + g/t, Maxwell (46) 
a = E€ + re, Kelvin - Voigt (47) 
a + Tid = E{e + T2e), 3 — parametermodel (48) 

with appropriate boundary conditions. All these models result in so-called 
creep and relaxation functions after the application of a unit stress or strain, 
as can be seen by the presence of the appropriate relaxation times r, ri and T2 
in the equations. Qualitative simulations of Maxwellian models have witnessed 
as one could expect increased tolerance to damage and crack meandering. This 
scenario arises due to the blunting of microcracks, and the fact that the local 
stresses are dissipated in the course of time. However, these models face certain 
challenges in dealing with boundary conditions, and stress- waves that originate 
from microfailures. 



4.5 Atomistic simulations 

In all materials the ultimate process leading to fracture is of course the break- 
ing of bonds at the level of individual atoms. This is true regardless of whether 
one considers metals with a crystal lattice, amorphous materials or fiber- 
matrix composites. In the fiber-matrix composites, it may be that the im- 
portant phenomena take place on a mesoscopic scale, however for crystalline 
solids (with dislocations present) it is clear that the fracture energy and the 
favored mechanisms of crack advancement depend in a fundamental way on 
atomistic effects. Since fracture is also crucially dependent on long-range elas- 
tic fields, one needs to consider processes on large, coarse-grained scales. The 
main question is thus how to bridge the length-scales from atomistic to contin- 
uum elasticity. A current trend is to consider multi-scale modeling in which one 
uses microscopic, even quantum-mechanical input to set up final simulations 
using the finite element method (see for example, Ref. [240]). 

A more fundamental problem is the following: since damage and deforma- 
tion on the atomistic level typically involve thermally activated or Arrhenius- 
processes, it is technically difficult to simulate a sample long enough to reach 
realistic conditions. Numerically accessible timescales are in fact in the range 
of picoseconds, much smaller than the normal laboratory scales on which frac- 
ture operates. One may of course try to circumvent this problem by using 
idealized systems like Lennard- Jones crystals [241-243]. This has been some- 
what popular for the last 20 years, and the crux of such simulations is that a 
prepared LJ crystal is subjected to shear or tension, leading finally to failure 
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Figure 29. Void initiation and subsequent growth in a binary LJ mixture. Each of the three rows 
represents five cross-cuts across a 3D system, while the rows themselves correspond to snapshots 
taken beyond the maximum stress (from Ref. [245]). 



exactly as in ordinary lattice models. In this case, however, one may incor- 
porate a finite temperature and rudimentary dislocations and their dynamics. 
Very large scale simulations by Abraham et al. [244] have demonstrated many 
fundamental phenomena in single crystal failure, namely, dislocation emission, 
sound waves and crack acceleration. As an example of recent application of 
this approach to phenomenological studies [245,246], we show in Fig. 1291 the 
formation of voids in a three-dimensional study of a binary amorphous LJ 
system. 

More elaborate potentials are then needed to make prediction on any real 
material, but the scenario gets even more complicated. For the elastic modulus, 
in many cases the available semi-empirical (to say nothing about full quantum 
mechanical ab-initio calculations) potentials allow relatively good quantitative 
accuracy. It is not clear, however, to determine the best possible way to define 
the condition for the failure of an atomic "bond" when a particular N-body 
potential is used. On a quantum mechanical level this concept makes no sense, 
and on a classical level one has to resort to ad hoc techniques. In addition, 
anharmonicity plays a significant role in such calculations. 

In spite of these difficulties, in recent times, there have been significant 
developments in using such atomistic simulations to study void growth in 
dynamic fracture of ductile metals, void growth ahead of the crack tip in silicon 
oxide and carbide, and crack propagation in single-crystal silicon [240, 247— 
250]. Another theme where progress is being made is in the context of failure 
of nanotubes [251]. However, such atomistic calculations are not yet at the level 
where they can be connected with mesoscale disorder, not even in the case of 
carbon nanotubes where fundamental issues still remain. The exception might 
be simulations of dynamic fracture, where in spite of the simplifications of 
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the potentials used one can nevertheless demonstrate fundamental phenomena 
[250]. We would however like to underline that even in such cases there have 
not been many studies of fracture as a stochastic phenomenon, with all the 
phenomenology that it involves. 

5 Statistical theories for fracture models 

The lattice models for fracture introduced in the previous section are rem- 
iniscent of other models studied in non-equilibrium statistical mechanics. It 
is thus tempting to apply standard theoretical tools to understand their be- 
havior and explore the analogies with other systems. This program turns out 
to be extremely complicated, due to a complex interplay between long range 
interactions, irreversibility and disorder. We discuss here the solution of sim- 
plified models (fiber bundles) and draw analogies between fracture and phase 
transitions. The fracture process is conveniently analyzed in the framework of 
interface depinning when disorder is weak, while the case of strong disorder is 
better viewed in the context of percolation. While these theoretical approaches 
provide a useful guidance, a complete understanding of the problem has not 
yet been reached. 

5.1 Fiber bundle models 

One of the simplest approaches to analyze fracture in disordered media is 
represented by the study of fiber bundle models [13,14]. These models were 
originally introduced to describe fibrous materials, schematizing the sample 
as a set of brittle fibers loaded in parallel. As in other lattice models, the 
failure threshold of each fiber is randomly chosen from a distribution. Next, 
one has to impose a rule for load redistribution after each failure. The simplest 
possibility is the case of an equal load sharing (ELS), in which each intact fiber 
carries the same fraction of the load. This case represents a sort of mean-field 
approximation and allows for a complete analytic treatment [14, 252-255] . At 
the other extreme lies the local load sharing model (LLS) where the load of a 
failed fiber is redistributed to the intact neighboring fibers [217,220,256-264] 
These simplified models serve as a basis for more realistic damage models 
such as the micromechanical models of fiber reinforced composites which take 
into account stress localization [256,258,259], the effect of matrix material 
between fibers [217,220,256,258,260-262], and possible non-linear behavior of 
fibers [265] . Other generalizations of the fiber bundle model include viscoelastic 
couplings [266-268], continuous damage [269], plasticity [270] and thermally 
activated fracture [271-273]. Fiber bundle models have been used in the past to 
address the macroscopic constitutive behavior, the reliability and size scaling 
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of material strength, and the avalanches of fiber breaks preceding ultimate 
failure [252-255]. 



5.1.1 Equal load sharing fiber bundle models. We consider first the case 
of ELS fiber bundles, in which N fibers of unitary Young modulus E = 1 are 
subject to an uniaxial load F. Each fiber i obeys linear elastic equation up to 
a critical load Xj, which is randomly distributed according to a distribution 
p{x). When the load on a fiber exceeds Xi, the fiber is removed. Due to the 
ELS rule, when n fibers are present each of them carries a load Fi = F/n and 
consequently a strain e = F/n. The constitutive law for ELS fiber bundles can 
be easily be obtained from a self-consistent argument. At a given load F, the 
number of intact fibers is given by 



n = N [1- / p{x)dx . (49) 




Rewriting Eq. (|49|) as a function of the strain, we obtain the constitutive law 

^ = 6(1-P(e)), (50) 

where P{x) is the distribution obtained from p{x). As a simple illustration, 
we consider a uniform distribution in [0, 1], so that Eq. (|5n)) becomes 

f = e{l-e), (51) 

where / = F/N . Similarly we can obtain the fraction of intact fibers p = n/N 
from Eq. (|49|) p = 1 — f / p which can be solved to yield 



p = ii + yr^)/2. (52) 

This equation shows that as the load is increased p decreases up to /c = 1/4 
at which p = pc = 1/2. For larger loads Eq. (|52|1 displays no real solution, 
indicating the onset of catastrophic failure. It is interesting to rewrite Eq. (|52|) 

as 

p = p, + A{f,-f)'/^, (53) 

with Pc = 1/2, fc = 1/4 and A = 1. In fact, this form is generically valid for 
most distributions p{x). To see this we rewrite Eq. H49() as / = x{l — P{x)) 
where x = f / p \s the load per fiber. Failure corresponds to the maximum Xc 
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of the left-hand side, after that there is no solution for Expanding close 

to the maximum we obtain / ~ + B{x — Xc)^, which then leads to Eq. (|53|1 
as long as the distribution is sufficiently regular (for a discussion of the limits 
of validity see Ref. [274]). This law also implies that the average rate of bond 
failures increases very rapidly before fracture: 

(54) 

The preceding discussion focused only on average quantities, but the pres- 
ence of a random distribution of failure thresholds necessarily leads to fluc- 
tuations. Due to the ELS rule, however, strength fluctuations are not partic- 
ularly interesting and vanish in the limit of large bundles. In particular, for 
any threshold distribution such that 1 — p{x) goes to zero faster than 1/x for 
X oo, the strength distribution is Gaussian with average fc = Xc{l — pixc)) 
and standard deviation a = Xcp{xc){l — p{xc))/^/N [14]. More interesting 
fluctuations are found in the precursors: the load redistribution after a fail- 
ure event can lead to avalanches of subsequent failures. This process can be 
treated analytically considering the sequence of external loads {-Pfc} at which 
some fiber fails [252-254]. After k — 1 failure event, the next fiber will break 
at a load = {N — k + l)xk, where the thresholds {x^} have been ordered. 
On average we have 

{Fk+i - Fk) = (1 - P{xk))pixk), (55) 

and close to global failure it vanishes as (Fk+i — Fk) ~ Xc — Xk- In addition, one 
can see that {{Fk+i — Fk) ) goes to a constant Hence, the external 

load {-Ffc} performs a biased random walk, and the bias vanishes as X ^ Xc 
When the load is fixed, an avalanche is defined by the return of the variable 
Fk to its initial value. For instance, given that the first fiber fails at load 
Fi, the avalanche proceeds until F^ > Fi. The process is then repeated from 
Fk, obtaining a series of avalanches whose size s is given by the first return 
time of the {F^} stochastic process. Thus the avalanche size distribution is 
equivalent to the first return time distribution of a biased random walk [254]. 
For a unbiased random walk, this distribution is given by D{s) = s~'^/^ so 
that the avalanches just before global failure decay as a power law. When the 
load is smaller than the critical load, the power law is exponentially cut off 
and the distribution becomes 

D{sJ)r. s-^/^eM-s/so), (56) 
where sq ~ (xc — x)^ ~ (fc — f)~^, diverging as f ^ fc is the limit N oo. 
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Eq. (|56|) refers to the avalanches at a fixed value of the load /. If instead 
we consider the integrated avalanche distribution obtained when the load is 
increased from zero to complete failure, we obtain 

Dint(s) ~ (57) 

as can be seen integrating Eq. (|56j) over /. Recently the difference between 
the integrated distribution, decaying as s~^/^, and the distribution sampled 
in a small bin, decaying as s"^/^, has been rediscovered and reinterpreted as 
a crossover indicating imminent failure [275]. In fact there is no crossover: the 
distribution at each load f < fc has an exponent 3/2 as stated Eq. 1)56^. but 
the cutoff prevents to see in practice the power law unless / ~ fc- 

5.1.2 Local load sharing fiber bundle models. The case of LLS represents 
another extreme case in which the effect of long-range stress field is completely 
neglected but stress enhancement around the crack is treated in the simplest 
way. Consider for instance a one dimensional series of fibers loaded in parallel 
with random breaking threshold from a distribution p{x). When the load on a 
fiber exceeds the threshold its load is redistributed to the neighboring intact 
fibers. Thus the load on a fiber is given by fi = /(I + k/2), where k is the 
number of failed fibers that are nearest neighbors of the fiber i and / = F/N is 
the external load [256]. Even for this apparently simple one dimensional model 
a closed form solution is not available, but several results are known from 
numerical simulations, exact enumeration methods or approximate analytical 
calculations. The analysis of the LLS model is extremely complicated and we 
thus list here the main results obtained, referring the reader to the relevant 
literature for the details [217,220,256-264]. 

Contrary to the ELS model, LLS fiber bundles normally exhibit non-trivial 
size effects as could be anticipated from general consideration of extreme value 
statistics. In particular, the average bundle strength decreases as the bundle 
size grows as 

/,~l/log(iV), (58) 

so that an infinitely large bundle has zero strength. The difference in the size 
effects between LLS and ELS is apparent in Fig. 1301 where the constitutive 
law of ELS fiber bundles is compared with that of an equivalent LLS model 
for different fiber bundle system sizes. The LLS follows the ELS law at small 
strains but fails before the ELS critical load. In addition, the failure point in 
the LLS model decreases as the lattice size is increased. In the limit of large 
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N, it has been shown that the strength distribution should follow the form 

W^(/) = l-[l-C(/)r, (59) 

where C(/) is a characteristic function, close to the Weibull form, but diffi- 
cult to determine exactly [256]. The existence of a limit distribution has been 
recently proved under very generic conditions for the disorder distribution in 
Rcf. [276]. This result implies that in LLS bundles there is no disorder induced 
transition from brittle failure, ruled by extreme value statistics, to a more grad- 
ual tough regime where fibers break gradually. The physics and mathematics 
of this problem is ruled by the determination of critical defects, and it is often 
so that one needs very large samples to meet some of these important critical 
defects (see Ref. [276] for a discussion about various LLS models). 

The avalanche behavior has also been studied both numerically [253] and by 
analytical methods [255]. While initially it was believed that LLS would also 
exhibit a power law avalanche distribution similar to ELS, it was later shown 
that the power law is only apparent and restricted to a small region. In fact, 
the integrated avalanche distribution is well approximated by 

Ant(s) ~ exp(-s/so) (60) 
where sq is independent on N. 

5.1.3 Generalizations of fiber bundle models. Fiber bundle models with 

ELS and LLS represent two extreme idealizations of fibrous materials that 
provide some insight on the general failure properties of disordered systems. 
These models have been improved in various ways, to obtain a more realistic 
representation of fibrous composites through a more detailed description of 
their mechanics and by including the effects of disorder on more complicated 
constitutive laws. There is a vast literature on the effect of the fiber arrange- 
ment, for instance two dimensional arrays rather than a one dimensional chains 
have been considered. Furthermore, it is possible to interpolate between LLS 
and ELS behavior through a long-range load transfer rule [269,277], or by 
a mixed mode in which a fraction of the load is transfered locally and the 
rest globally [278]. In both cases one observes a crossover between local and 
global behavior as a function of the parameters of the model. In particular, in 
Ref. [277] the authors use a load transfer law decaying where r is the 

distance from the broken fiber, and study the failure process as a function of 7. 
For 7 close to zero, they recover the ELS result, while the LLS results is found 
for large 7. The two regimes are separated by the point 7c — 2. This model 
is interesting because it makes contact with lattice models for fracture where 



69 




stress is generally transferee! in a non-local, but geometrically dependent, fash- 
ion. Interestingly, for the two dimensional random fuse model the current is 
redistributed with a law decaying as exactly the crossover point between 
local and global behavior. We will elaborate further on this relation in the 
next section. 

The variable load transfer rule also allows to study the patterns of damage 
using 7 as a control parameter, as does the version where perfect plasticity 
is included [270]. In the latter case, with LLS one can have a growing crack 
which turns out to be compact, with a rough perimeter. The connected cluster 
of broken fibers can also percolate forming a particular case of a percolation 
transition. This is slightly different from the damage patterns observed in more 
realistic - even RFM - models, which are isotropic. Last, in the partly plastic 
FBM the avalanche size distributions of bursts deviate from their LLS/ELS 
behaviors [270]: the scale-free nature vanishes with ELS, while there is a special 
value of remnant load in a yielded fiber, that seems to induce a power-law even 
with LLS. 

Next, we discuss in more detail some time dependent generalizations of fiber 
bundle models. In viscoelastic fiber bundles [266-268] each fiber obeys a time 



70 



dependent constitutive equation of the type 

f = Pe + e, (61) 

where / is the external load, /3 is a damping coefficient and the Young's mod- 
ulus is unitary. In a creep test (i.e. constant load test), one thus finds that the 
strain evolves according to 

e(t) = /(I - exp(-V/3)) + eo exp(-V/3) (62) 

where €q is the initial strain. Fibers fail when the local strain e exceeds a ran- 
dom threshold Xi, taken from a distribution p{x), and its load is redistributed. 
Here we discuss the ELS case, the generalization to LLS being straightforward. 
The asymptotic behavior of the model is the same as in the time independent 
model since the Hooke's law is eventually valid when t S> /3, but the order 
of individual failures is different. The time dependent constitutive behavior of 
the model is obtained by solving the self-consistent equation 

/ = (/3e + 6)(l-P(e)), (63) 

reducing to Eq. H50|) in the steady state regime. A typical example of the time 
evolution for different values of the loads is reported in Fig. 1311 From this we 
can expect that the bundle will eventually fail for f > fc and resist otherwise. 
The time to failure can be evaluated from Eq. (|63() and is found to diverge as 

~ (/ — /c)~^/^ as / — > fc- Notice that in the LLS version of the model one 
still observes a transition but t f does not diverge as / ^ [267] . 

Another interesting time dependent generalization of the fiber bundle model 
considers thermally activated fracture [271-273]. The additional ingredient 
that is added to the model, typically in the ELS approximation, is the presence 
of an uncorrelated Gaussian noise term 7/(t) acting independently on each fiber. 
In particular the noise distribution is given by 

pW = ^exp(-|i), (64) 

where T is the temperature in appropriate units. Contrary to the viscoelastic 
case, the noise has the effect to induce fracture even in subcritical conditions, 
when / < /c. It is instructive to consider first a bundle without disorder, in 
which p{x) = 5{x — 1). Even in this simple case, it is not possible to obtain 
a general closed form for the failure time, but in the low temperature limit 
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Figure 31. The time evolution of the strain of a viscoelastic fiber bundle for different values of the 
load /. For f < Fc the strain accelerates up to a time to failure tf (from Ref. [268]). 

T « 1 — / case, one can obtain an approximate expression 

In the presence of disorder the problem becomes clearly more complicated, 
but there is an indication that in the low temperature case, the failure time 
follows a law similar to Eq. (|65|) with an important difference: the temperature 
T is replaced by an effective temperature T^jf = T + 0, where depends on 
disorder. Hence the net effect of disorder is to enhance thermal activation. It is 
remarkable that an analogous result was found in Ref. [23] for crack nucleation 
in a disordered medium. 



5.2 Statistical mechanics of cracks: fracture as a phase transition 

The idea that there is a relation between fracture and phase transitions has 
a long history. For instance, the Griffith theory of fracture is very similar in 
spirit to the classical theory of nucleation in first order phase transitions. In 
bubble nucleation, a critical droplet will form when the loss in free energy due 
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to the bulk forces exceeds the increase in the interfacial energy. Similarly frac- 
ture occurs if the external stress prevails over the resistance at surface of the 
crack. Further analogies come from the scaling behavior observed in fracture 
experiments, such as for the crack roughness and the acoustic emission distri- 
butions. In the following we provide a basic introduction to phase transitions 
and discuss its relation to fracture using lattice models as an illustration. 



5.2.1 Generalities on phase transitions. Phase transitions are character- 
ized by changes in the internal symmetries of a material as external control 
parameters are varied. Familiar examples are the melting of a crystal, or the 
ferromagnetic transition in a magnet. In the first example, we have an abrupt 
first-order phase transition, with latent heat, coexistence and no precursors, 
while the latter is a continuous second-order transition. Notice that one can 
also observe a first or a second order transition in the very same system, de- 
pending on the external conditions, the typical example being the gas-liquid 
transition that is first order at low temperatures and pressures and becomes 
second order in a particular point of the PT diagram. 

To be more quantitative, one typically distinguishes a phase by an order pa- 
rameter whose value is related to the internal symmetries of the system. For 
instance, a ferromagnet acquires a non-zero magnetization when the spin ro- 
tational symmetry of the paramagnetic phase is broken. Since the transition is 
continuous, the magnetization vanishes as the transition is approached. This 
is unlike first order transitions, where the order parameter is discontinuous 
across the transition. This case is typically associated also with metastability: 
when the transition point is reached, the system needs some time to cross 
into the new equilibrium state. It does so by nucleation: small regions of the 
new phase are formed, until they become large enough to transform the sys- 
tem completely in a rapid event. In second order transitions, the buildup of 
the new phase instead occurs before the transition point is reached through 
the formation of larger and larger correlations that span the entire system 
exactly at the critical point. This process is encoded in scale invariance: ther- 
modynamic quantities around the critical point are described by homogeneous 
scaling functions. 

To be more concrete and set up the notation, we consider the example of 
a uniaxial ferromagnetic system. The control parameters, in this case, are 
the temperature T, the magnetic field H and the average magnetization M, 
which is the order parameter of the system. At high temperatures, fluctuations 
dominate and the magnetization is zero on an average. As the temperature is 
lowered, the local magnetization becomes more and more correlated until the 
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critical point T = Tc is reached. The correlation length ^ then diverges as 

(66) 

where v is a critical exponent. The average susceptibility x = dM/dH also 
diverges, indicating that the system is extremely prone to changes, and it does 
so as 

X-^{T-T,)-\ (67) 

The order parameter M becomes non-zero in the ferromagnetic phase and 
close to the critical point follows 

M ~ (T - Tc)^, (68) 

defining another critical exponent. The presence of a magnetic field destroys 
the transition, forcing a particular direction for the magnetization, and thus 
defines another scaling law at T = Tc 

(t> ~ H^'^. (69) 

We have defined a certain number of critical exponents, which are however not 
all independent, as can be shown by a more detailed analysis. Here, we do not 
want to go further into these details but just highlight the fact that second- 
order phase transitions are associated with scaling and a diverging correlation 
length ^. 

A ferromagnet can also be used to illustrate the idea of nucleation, which is at 
the core of first order transitions. Imagine that the temperature is fixed below 
the critical point and we instead vary the field H. When the field is suddenly 
increased from a negative to a positive value the magnetization should also 
change sign, but if the value of the positive field is not so strong then the 
state with negative magnetization is still metastable. Thus small local changes 
of the magnetization will be reabsorbed, while large sudden changes will be 
less favorable for entropic reasons. The system will spend a considerable time 
in this metastable state before a large enough change occurs. To be more 
precise, the free energy cost of a spherical droplet of size I is the sum of two 
contributions: a bulk decrease due to the magnetic field and a surface energy 
increase: i.e. 

AF = -Att^HI^ + ttJI^, (70) 
where J is the surface energy, li I > l^. the droplet is unstable and a further 
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increase in droplet size would result in a decrease in its free energy. Thus 
the nucleation time can be estimated as the time needed to form the first 
unstable droplet. The analogy with Griffith theory is evident and one could 
thus associate crack growth with a first order transition from an elastic body 
to a fractured solid driven by stress. 

The difference between first and second order behavior becomes more nu- 
anced when nucleation takes place close to a spinodal point, corresponding to 
the field at which the metastable state becomes unstable and thus the sys- 
tem is forced to reverse. Strictly speaking, the spinodal point only exists in 
mean-field theory, but its signatures can be observed even when interaction 
is long-ranged. The theoretical description of homogeneous spinodal nucle- 
ation is based on the Landau-Ginzburg free energy of a spin system in the 
presence of an external magnetic field. When the temperature is below the 
critical value, the free energy has the typical two- well structure as depicted 
in Fig. 1321 In the presence of an external magnetic field, one of the wells is 
depressed with respect to the other, which represents therefore the metastable 
state. The system must cross a free energy barrier to relax into the stable 
phase. When the external field is increased, this nucleation barrier decreases, 
eventually vanishing at the spinodal. The approach to the spinodal is charac- 
terized by scaling laws, analogous to critical phenomena, but the reversal is 
discontinuous and the transition is first-order. The magnetization or the order 
parameter M scales with the external field H as 

M - Ms-- [Hs- Hf/"^ (71) 

where Mg and Hg are the values of the order parameter and the field at the 
spinodal. This law implies a divergence of the quasi-static susceptibility 

X ^ ^ ~ (i^s - H)--^- 7 = 1/2 (72) 

The fluctuations in the order parameter can be related to suitably defined 
droplets, whose sizes turn out to be power law distributed with an exponent 
r = 3/2. For finite-dimensional short-range models, this mean-field picture is 
expected to fail, since the system will nucleate before reaching the spinodal 
point. Thus no scaling would be observed in this expected for first- 

order transitions. 

5.2.2 Disorder induced non- equilibrium phase transitions. We have dis- 
cussed above the properties of nucleation and phase-transitions in thermally 
activated homogeneous systems. In the case of fracture, however, we often deal 
with a mechanically driven disordered system. An interesting analogy can be 
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made with the zero-temperature random-field Ising model (RFIM) proposed 
by Sethna et al. [279-281] in the context of magnetic hysteresis. In this model, 
a set spins Sj = ±1 are assigned to the sites z of a lattice. The spins interact 
with their nearest-neighbors by a ferromagnetic coupling J and are subject 
to an external field H. In addition, to each site of the lattice is associated a 
random field hi taken from a Gaussian probability distribution with variance 
i?, such that p{h) = exp(— /t^/2i?^)/-\/27ri?. The Hamiltonian reads 

E = -Y,JsiSj-Y,{H + hi)si, (73) 

where the first sum is restricted to nearest-neighbors pairs. 

The system is started from a saturated negative state and the external field 
is increased from a large negative value. At each time the spins Sj take the 
sign of the local field 

hf^'^-'^^jY.^.+h + H, (74) 
' i 

where the sum is over the nearest neighbor sites j. When a spin i flips, the 
effective field of the neighbors j is increased by 2J. This can lead the effec- 
tive field to change sign inducing an avalanche. For small values of the order 
parameter A = R/J, flipping a few spins generates a big avalanche whose 
size is comparable to the system size, leading to a discontinuous magnetiza- 
tion reversal. On the other hand, a large disorder prevents the formation of 
large avalanches and the magnetization reversal is a smooth process. The two 
regimes are separated by a critical point Ac, where the avalanches are dis- 
tributed as a power law. The behavior of the model in the (A, H) parameter 
space is very similar to that of standard equilibrium phase transitions as it is 
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shown in Fig. 1331 In particular, all the usual exponents discussed in Sec. 15.2.11 
can be defined replacing T with A. 

In the mean-field theory, the approach to the instability He, is characterized 
by the same scaling laws and exponents of spinodal nucleation, as reported in 
Eqs. (|71() ~()72 |) [280]. In addition, the avalanche size distribution is described 
by a scaling form 



with r = 3/2 and k = 1. It is worth noting that these scaling exponents co- 
incide with the mean-field exponents for the distribution of droplets in homo- 
geneous spinodal nucleation. From these studies it appears that the behavior 
of thermally activated homogeneous spinodal nucleation is similar to the ap- 
proach to the instability in disordered systems driven at zero temperature. 
However, one should bear in mind that for a given realization of the disorder 
the dynamics is completely deterministic in the second case. Concepts such as 
metastability and nucleating droplets are formally not defined in this context. 
This picture is only valid in mean-field theory or when interactions are long- 
ranged. As in conventional equilibrium phase transitions, when interactions 
are short-ranged nucleation takes place before the spinodal line and scaling 
laws can be observed only at the critical point. 

5.2.3 Phase transitions in fracture models. The interpretation of fracture 
as a phase transition is a highly controversial question. While at first glance 
the analogies are striking, a precise relation is hampered by several difficulties 
and no consensus has been reached in the literature on this subject. A series 
of studies consider the case of crack nucleation in homogeneous media, elabo- 
rating on the analogy between the classical theory of nucleation and Griffith 
theory. Selinger et al. [241-243] have used mean-field theory and numerical 
simulations to show that a solid under stress is in a metastable state, becom- 
ing unstable when the external stress reaches a spinodal point. This is best 
illustrated with a simple linear chain of particles. It breaks down once one of 
the pairs gets separated due to a thermal fluctuation that makes the inter-bond 
to fail. Rundle and Klein [282] have proposed an equation for the growth of a 
single crack, deriving scaling laws for spinodal nucleation. The nature of the 
nucleation process in a stressed solid was further analyzed in Refs. [283, 284] 
using Monte Carlo simulations. More recently, it was shown more rigorously 
in the framework of elastic theory, that the point of zero external stress cor- 
responds to the condensation point in gas-liquid first order transitions [285]. 
It is important to remark that these studies deal with the thermally activated 
fracture where quenched disorder is not considered. The presence of disor- 
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Figure 33. The phase diagram of the non-equiUbrium RFIM. The dashed Hne represents the 
location of the spinodal, which is obtained from the mean-field theory. When the interactions are 
short-ranged the reversal occurs on the solid line through a first-order transition. The first-order 
line ends at a critical point, above which the reversal is continuous. 



der, in the form of vacancies or microcracks, strongly affects the nucleation 
process [242,283,284], providing in some cases an effective disorder-induced 
correction to the temperature [23] as pointed out above in Sec. 15.1.31 

Several authors have tried to relate the scaling laws observed in fracture mor- 
phologies and acoustic emission distributions to an underlying critical point, 
basing their analysis on experiments [121,124] and models [253,286-290]. A 
natural starting point is provided by a lattice model, which could be mapped 
into other models, typically drawn from ferromagnetism, where the presence 
of a phase transition is well established [286,287,289,290]. It is instructive to 
consider first the fiber bundle model, since exact results are available and the 
nature of the hypothetical phase transition can be clearly identified [286,287]. 
To make contact with the RFIM, we can assign to each fiber a spin variable 
Si, whose value depends on whether the fiber is intact (s = —1) or broken 
{s = 1). The effective field is formally given by 



hf^^=F,{{sj},F)-Xi 



(76) 
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where Xi is the fiber threshold, Fi is the load on the fiber i, which depends on 
the external load F and on the state of the other fibers. The dynamics is now 
similar to the one of the RFIM since each "spin" follows the local field (i.e. 
it breaks when Fi > Xi). Thus changing the disorder strength and the nature 
of the load sharing we could obtain a phase diagram of the type depicted in 
Fig. 031 

In the case of ELS fiber bundles, the load on each fiber is given by Fi = 
2F I (N — Y2j while for LLS Fi is a complicated function of the state of the 
neighboring fibers: if a fiber is close to a crack of length k it carries a load of 
Fi = F{1 + k/2). The crucial difference between these load sharing rules and 
ferromagnetic interactions is that for fibers the increase of the effective field 
due to other fiber breaking grows as the fracture process progresses, while in 
the RFIM this quantity is always given by 2 J. For instance in the case of ELS, 
the breaking of a single fiber increases the effective field of the others by a 
quantity AF = F/{n{n — 1)) which increases as the fraction of intact fibers n 
decreases. Similarly in the case of LLS each time we break a fiber, we increase 
the effective field of the neighbor by a quantity AF = kF/2 which becomes 
larger as more and more fibers break. 

This difference has a profound effect on the phase diagram, because it im- 
plies that interactions always prevail over disorder. Thus, as the system size 
increases, the system will necessarily hit the low disorder first-order transition 
line. The only exception is provided by the extreme case of a non-normalizable 
disorder distribution, where we can always find some strong enough bond to 
resist the load amplification [274,288]. In more general cases, however, it is 
not possible to fine tune the ratio between disorder width and interactions as 
for the RFIM in order to find a critical point where the two effects are per- 
fectly balanced. Thus the "phase transition" observed in fiber bundle models 
is always first-order, with a catastrophic nucleation event preceded by some 
small precursors. The scaling observed in the ELS models is the one associated 
with a spinodal instability, similar to what is found in the low disorder phase 
of the RFIM in mean-field theory. When interactions are local, such as for the 
LLS model, the spinodal line can not be reached and the transition is simply 
first order. With this idea in mind we can reinterpret the results presented in 
Fig. 1301 as the range of interaction is increased nucleation occurs closer and 
closer to the spinodal point, as in conventional first-order phase transitions. 

Having settled the question of fiber bundle models, we can now turn our 
attention to more complicated lattice models, such as the random fuse model. 
The attempts made in the past to map fracture models into spin models are 
typically based on the lattice Green function formalism: removing a bond i in 
the lattice implies that the load in the other bonds has to be changed by a 
quantity Gij (the Green function). Thus we can rewrite the effective field in 
Eq. (|76j) in terms of Gij and obtain a suitable spin model. For instance, in the 
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random-fuse model Gij can be approximated as the dipolar kernel, but only 
in the limit where only a few bonds are broken. It is interesting to remark 
that in these conditions the Green function decays as r^^ in two dimensions, 
which is exactly the crossover point between local and global behavior in the 
FBM. This could imply that the spinodal point, inaccessible under LLS, could 
be reached in the RFM, as is also suggested by the simulations discussed in 
Sec. EH 

One should take these considerations with due care, since for a generic dam- 
age configuration, finding Gij amounts to solving the full problem (i.e. the 
Kirchhoff equations for the random fuse model) and is thus practically im- 
possible. Thus we can rigorously map fracture models into spin models only 
in the dilute limit, which is not, however, the limit of interest in fracture. 
Extrapolating the results of spin models outside their regime of validity leads 
to inconsistent results, because the stress amplification at the crack tips is 
not properly captured by the dilute Green function. For instance, the authors 
of Ref. [289] treat in mean-field theory the spin model obtained from Green 
function formalism. Neglecting the load amplification in their analysis (i.e. the 
load transfer is chosen to be independent of the state), they obtain a spurious 
phase diagram similar to that obtained for the RFIM. Similar problems affect 
the analysis presented in Ref. [290]. 



5.3 Crack depinning 

We have discussed in Section IHl that experiments in several materials under 
different loading conditions, lead to a rough crack surface described by self- 
affine scaling. The simplest theoretical approach to the problem identifies the 
crack front with a deformable line pushed by the external stress through a 
random toughness landscape and the trail of the line yields the fracture surface 
[91,97]. The roughness of the line comes from the competition between the 
deformations induced on the line by the materials inhomogeneities and the 
elastic self-stress field that tries to keep the front straight. This problem has 
been studied in full generality and the basic phenomenology applies to a variety 
of systems including domain walls in ferromagnetic and ferroelectric materials, 
dislocations gliding through solute atoms, flux lines in superconductors and 
others. 

For simplicity, it is instructive to consider first the case of a planar crack 
front, corresponding to the experiments reported in Refs. [83-85]. In this case 
we can schematize the crack as a line moving on the xy plane with coordinates 
{x, h{x, t)) (see Fig. lSl)) [291-293]. An equation of motion for the deformed line 
position is obtained by computing, from the theory of elasticity, the variations 
to the stress intensity factor induced by the deformation of the front. In the 
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quasistatic scalar approximation, this is given by [294] 

Km^m =K,j dx'^^^^^^^^^^, (77) 

where Kq is the stress intensity factor for a straight crack. The crack deforms 
because of the inhomogeneities present in the materials, and these inhomo- 
geneities give rise to fluctuations in the local toughness Kc{x, h{x,t)). These 
ingredients can be joined together into an equation of motion of the type 

dh 

r-^ = Ke^t + K{{hix, t)}) + Kcix, h{x, t)), (78) 

where F is a damping term and Kgxt is the stress intensity factor correspond- 
ing to the externally applied stress [291-293]. Eq. (|78|) belongs to a general 
class of interface models which have been extensively studied numerically and 
analytically. For low stress the crack is pinned, and there is a critical threshold 
Kc above which the crack advances at constant velocity 

V ~ (EText - K,f, (79) 

where /? is a scaling exponent, whose value is estimated numerically to be 
P ~ 0.68, while a renormalization group analysis yields /3 = 2/3 [295]. Close 
to the depinning transition crack motion becomes correlated, with avalanches 
of activity. The distribution of avalanche lengths I follows a power law distri- 
bution with a cutoff 

p(/)~r''/(//e), (80) 

where k ~ 1 and the correlation length diverges at the transition as 

e ~ {Kext - K,)-\ (81) 

with an exponent estimated as ~ 1.52 [293]. In addition, the crack front is 
self-affine, with a roughness exponent given by 

C = 1 - (82) 

Numerical estimates of the roughness exponents range between = 0.34 [291, 
292] to C = 0.39 [296]. 

The results of the planar crack propagation model can be compared with the 
experiments reported in Refs. [83-85]. It is apparent that there is significant 
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discrepancy between theory and experiments, since the roughness exponent 
is considerably larger in the latter case. Several mechanisms have been pro- 
posed in the literature to account for this discrepancy, but none of them so 
far is entirely convincing or commonly accepted. The arguments involve cor- 
related disorder [297], elastodynamic effects [293], crack front waves [229] and 
microcrack nucleation ahead of the main crack [298-301] . 

The more general problem of three dimensional crack surface roughening 
can also be recast as a depinning problem, with similar difficulties [230]. It 
is possible to derive an equation of motion, based on linear elasticity, for the 
evolution of the crack surface. With respect to the planar crack case, we have 
an additional equation for the out of plane component h±{x,y) of the crack, 
while the in-plane component h^\{x) follows Eq. (|78|). The equation for h{x,y) 
and the resulting roughness exponent depend on the particular fracture mode 
imposed in the experiment. Without entering into the details of the derivation, 
which can be found in Ref. [230], we report here the result. In mode I the out 
of plane roughness is only logarithmic, while an exponent (± = 1/2 is found 
in mode III. Both results are quite far from the experimental measurement 
C ~ 0.78, although the mode III result could be in line with the value C = 0.5, 
sometimes measured at small scale (consider also Ref. [92]). Unfortunately, 
experimental results are normally obtained under mode I loading conditions, 
at least effectively, while mode III cracks tend to be unstable. 

Given the results discussed above, we would conclude that crack rough- 
ening is probably more complicated than a simple depinning problem of a 
linear elastic line moving in a disordered landscape. Nonetheless the depin- 
ning picture is quite appealing, since it relates the observed scaling behavior 
to a non-equilibrium critical point. To obtain a quantitative explanation of the 
experiments, it is necessary to analyze in detail further modifications of the 
basic depinning problem, including for instance the combined effects of the 
anelastic fracture process zone, the formation of microcrack, plastic flow, or 
elastodynamics. Among these proposals the most promising pathway for future 
investigation is probably represented by the effect of non-linear elasticity. 

The depinning of a non-linear long-range elastic line has been recently stud- 
ied by renormalization group for contact line depinning, a problem that is for- 
mally equivalent to a planar crack [302] . It was shown that adding a non-local 
KPZ-type non-linearity ^ to Eq.[7Hl raises the roughness exponent to C — 0.45. 
This value is a little closer to the experimental result, but needs further confir- 
mation, since it is based on a one-loop e expansion. A more striking similarity 
is obtained if we consider a non-linear, but local, elastic kernel in Eq. 1781 In 
this case, numerical simulations indicate that the exponent is C = 0.63 [303] 



^The standard KPZ term is A(V/i)^, a non-local KPZ term is written in Fourier space as Xq°'\h{q)\^ 
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which is in perfect agreement with the experiments. This universahty class cor- 
responds to the pinning of the front by a directed percolation path [304,305]. 
In order to accept this explanation, we have to understand the mechanism 
responsible for the screening of long-range forces, resulting in a short-range, 
but non-linear, kernel. In this respect, recent experiments on amorphous Silica 
show that the roughness correlation length ^ is of order of the size of the FPZ 
(,FPZ [95,96,306]. By definition, linear elasticity is not applicable inside the 
FPZ, where the crack roughening process would takes place. Hence, this obser- 
vation could be used to justify why the effective kernel is local and non-linear, 
accounting for the observed value of (. While this appears to be a plausible 
explanation, further work is needed to substantiate this claim. 

5.4 Percolation and fracture 

5.4.1 Percolation scaling. Adding randomness to brittle lattice models of 
fracture can be done in many ways, the simplest being the random removal of 
a fraction of the elements. This process is best understood using the concept 
of percolation, a particular second-order phase transition that exhibits the 
full spectrum of usual characteristics, such as scaling and critical exponents as 
discussed in Sect. 15.2711 [307]. A percolation transition signifies that a structure 
becomes connected for certain values of the control parameter (i.e. the fraction 
of intact bonds). The usual example, demonstrated in Fig. I35( is that of bond 
percolation in a two-dimensional square lattice, where a fraction p of the bonds 
have been removed. The question is now whether the remaining ones form a 
connected (or spanning) cluster that extends through the lattice, connecting 
for instance the two horizontal edges. In this particular example, the spanning 
cluster is found, in the limit of large lattices, for p < pc = 1/2. As in other 
second-order phase-transitions, the location of the critical point (i.e. the value 
of pc) depends on the type of lattice, but the critical exponents are universal. 

An important quantity for fracture in the percolation context is the order 
parameter, usually defined as the probability 11 that a bond belongs to the 
spanning cluster. Below the percolation threshold pc, the order parameter 
scales as 



with /3 = 5/36infi = 2, while 11 = for p > pc- The essential part of 
describing percolation as a second-order phase transition is that there is a 
divergent correlation length ^ scaling as 
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Figure 34. The geometry of a deformed planar crack 



with V = 4/3 in (i = 2. This correlation lengtli can be measured analyzing 
the statistics of the clusters of connected bonds for y < pc- The cluster ra- 
dius distribution is a power law up to a cutoff lengthscale ^. For lengthscales 
larger than ^ (for the system size L S> ^ or for distances between two points 
r ^ S^) the system is out of the critical regime and is either connected or dis- 
connected. For scales such that r < ^, the physics of percolation is governed 
in general simply by the singly connected or "red" bonds. They are - in the 
bond percolation language - those that make the spanning cluster connected, 
so that removing any of them also destroys connectedness, or cuts the continu- 
ous, spanning path across the sample into two. The sensitivity to the lattice to 
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Figure 35. The backbone of a bond percolation problem in 2d: only the bonds belonging to it are 
depicted. The "red bonds" arc marked, naturally, with red. So-called busbar boundary conditions 
are used, which makes the backbone more dense close to the boundaries. 

small perturbations, like the removal of a red bond, is a signature of criticality. 

Red bonds also play a crucial role on various important properties. For 
instance, if the bonds are taken to be fuses in a RFM, then red bonds will 
all carry the same current. Thus the conductivity (or equivalently the elastic 
stiffness in tensorial models) and strength are going to relate to the properties 
of such bonds. This allows to write down in the critical regions the scaling 
Ansatz for the conductivity S 

^^{p-PcY^fperciC/L'^n (85) 

such that for L ^ ^ or p large enough but still close to the critical regime, the 
scaling function / becomes a constant. We have formulated Eq. H85|) in terms 
of a (general) transport exponent tj], which can also describe any particular 
physical property from conductivity to maximum current to strength, if the 
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bonds are taken to be beams with certain elasticity properties. The ratio of the 
conductivity exponent ts to v is for instance known to be ts/z^ = 0.9826. . . 
[308] in 2d, whereas in 3d the ratio is about 2.3 [309]. These questions are 
often analyzed either numerically via transfer matrix calculations on strips, 
by renormalization group arguments, or by real-space renormalization, and it 
appears that the numerical value of most any exponent to is non-trivial in the 
sense that it is not simply related to u and the physical dimension. 



5.4.2 Variations of the percolation problem. Interesting physics is en- 
countered in several separate limits, which concern the nature of the inter- 
action Hamiltonian imposed on the bonds. The difference between clastic and 
scalar /electric percolation should also be visible in the strength case. At the 
minimal level, this follows from the scaling of the stiffness and conductivity 
with p — pc- A fundamental special case arises if the angular or bond-bending 
part is taken to be zero, exactly, so that only central forces remain. Then 
the main effect is that single-connectedness is no longer sufficient since simple 
Maxwell-counting demonstrates that the degrees of freedom are not balanced 
by the constraints at hand. Several often-used lattices become unstable (e.g. 
simply cubic ones, like the square one) against deformation. The critical be- 
havior is now called rigidity percolation, and its nature is still a not completely 
resolved question, e.g. with regards to its order in two dimensions - whether 
it is simply second-order similar to percolation or not [310]. For studies of 
central-force system fracture in diluted systems the implications are not quite 
clear either. With increasing damage, at some point, the system will have a 
vanishing (linear) elastic modulus, while still being geometrically connected. 

Another example of percolation-related processes in the fracture context is 
given by continuum percolation, where one considers the transport proper- 
ties of a medium with a set of voids or defects. The difference to the lattice 
case is now due to the fact that the "bottlenecks" of the spanning cluster 
can have a (naturally) varying distribution of properties. A typical example is 
the so-called Swiss Cheese -model, where one introduces spherical or circular 
defects [311,312]. From geometrical considerations, one can now compute the 
distribution of minimum widths for the narrowest ligaments that one finds 
on the spanning cluster. The result is that transport properties are changed. 
In this example, the conductivity is expected to follow S ~ (p — Pc)'^"" with 
Tsc = {d — 2)z/ -|- max[9, (1 — a)^^]. Here 6 is the resistance exponent of nor- 
mal percolation, describing what happens between two terminals on the same 
cluster. The exponent a depends on model and dimension, and describes the 
power-law distribution of the conductivities of the bottle-necks. One can de- 
rive the failure properties easily by noting that all the red bonds carry the 
same current, but the result depends on what one assumes about the failure 
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thresholds (see also Ref. [313]). 

Finally, an opposite limit to the diluted case is given by the infinite disorder 
case, in which the medium fails following two rules: one picks always the bond 
with the smallest threshold, and excludes in the process those that do not 
carry current. Evidently in the beginning this is just dilution, as in geometric 
percolation, but at some point the screening will start to play an increasing 
role [314] . In this context and in studying fracture close to Pc it might be helpful 
to use recent algorithmic developments in pruning percolation clusters [315]. 
These allow to consider only the backbone, and to identify the red, critical 
bonds, with polynomial time scaling in L and could thus be used to obtain 
starting configurations for simulations. 



5.4.3 Strength of diluted lattices. Above the critical region, when the lat- 
tice is connected, one can consider the scaling of strength as a function of L 
and p. For simplicity, we restrict our discussion to the scalar case. As noted 
above, the random removal of bonds naturally leads to a macroscopic reduction 
of the conductivity from the homogeneous system value. Figure OHl shows an 
example of the behavior in two dimensions, with fixed L [316]. The crucial ar- 
gument that relates the strength ac{p,L) to the extremal statistics arguments 
underlying the scaling of brittle fracture is due again to Duxbury, Beale, and 
Leath [316,317]. The most critical defect - this holds also for the elastic case - 
is a linear one, a row of n missing bonds in two dimensions. In higher dimen- 
sions, the corresponding situation is given by a penny-shaped crack and the 
fuse current will scale as iup ~ io{l + k^n^^^'^'^^^^). This results from noticing 
that the total current increase will scale as n^/^, but it is shared by a number 
of fuses that depends on dimension (e.g. d = 2 in 2d). 

Forgetting about rarer configurations where two such cracks are found next 
to each other, the probability of having a defect of given length is clearly 
proportional to P„ ~ P^(l — p)"', times a prefactor which is L-dependent from 
the row length. Thus a connection can be established non-rigorously to the 
basics of extremal statistics distributions outlined earlier: the critical defect is 
the biggest found in a sample of transverse size (L, usually) and consisting of 
(again) L independent subvolumes. 

From the damage mechanics viewpoint it is useful to take note that in this 
variant of the RFM the fracture process zone is practically non-existent; thus 
there is no corresponding lengthscale S,fpz 3> 1 that might depend on system 
size. Furthermore, for any p the damage or number of fuses that fails before the 
maximum current is very small. Thus in these systems there is no "nucleation" , 
the largest defect in the original system dictates the strength. 

In particular, the defect population with dilution or percolation disor- 
der is exponential, and thus the limiting distribution will be of the double- 



87 



it X 



1,* 



Figure 36. The scaling of both conductivity and strength in 2d RFM as a function of the bond 

probability p (from Ref. [316]). 

exponential or Gumbel-form, of extremal statistics [318]. One can write for 
the integrated distribution 



P{i) = 1 - exp (-cL<^exp [-A;^^ ^/"]), 



(86) 



where a is a constant that depends on p (recall again this is for percolative 
disorder), as are c and k. For the average strength it follows that there is a 
logarithmic size-effect, {ac{L)) ~ 1/logL" where a ~ 0{1) and is difficult to 
establish rigorously. 



5.4.4 Crack fronts and gradient percolation. An interesting variation of 
the percolation problem is gradient percolation in which the symmetry is bro- 
ken by letting the occupation probability of bonds to be a function of the 
vertical coordinate y so that p = p{y) [319]. In the simplest case, p{y) is a 
linearly decreasing function interpolating between p = 1 and p = 0, but other 
choices do not change the basic physics. This model has applications in several 
mean-field -like fracture models in which the gradient is imposed into damage 
by an external displacement or potential field [298-300]. 

The essential feature of gradient percolation is that an occupied region is 
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Figure 37. Example of a 2d fuse network for p = 0.75. What is the largest linear defect in such a 

system? 

separated from an empty region by a corrugated front. The front occurs around 
the region where p{x) ~ Pc and displays a self-similar geometry at short scales, 
crossing over to a flat profile at larger scale. The crossover scale is set by 
the value of concentration gradient V^* around Pc- Using scaling arguments it 
is possible to show that 



a 



(87) 



where is the correlation-length critical exponent of percolation defined in 
Eq. (|5Hl. The width of the separation front scales linearly with and thus 
increases as the gradient becomes weaker. For lengthscales smaller than the 
front is essentially equivalent to the perimeter of a percolation cluster and has 
thus a fractal dimension D = 4/3 (or D = 7/4 depending on the definition of 
the front) [320]. An example of a gradient percolation cluster is reported in 
Fig.M 

The relation between gradient percolation and fracture is based on the study 
of models of planar cracks. An approach alternative to the crack line depin- 
ning scenario discussed in Sec. 15.31 considers damage accumulation in a stress 
gradient. Thus the crack is not considered as line, but as the separation front 
between damaged and intact areas, allowing for microcrack formation and co- 
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alescence even ahead of the main crack. The question is whether this effect can 
account for the discrepancy between the experiments and the crack depinning 
predictions. 

Different types of lattice models have thus been simulated with boundary 
conditions inducing a stress gradient. Simulations of a two dimensional RFM in 
the planar crack geometry (i.e. considering two conducting plates connected by 
fuses) yield a gradient percolation front [298] . A beam model studied in a sim- 
ilar geometry, recovers instead the crack line depinning exponents, indicating 
that damage nucleation is irrelevant for crack roughness in that model [299]. 
On the other hand, a long-range scalar model, simulating planar cracks in 
a three dimensional geometry, yields a larger exponent ~ 0.6 which could 
explain the experiments [300,321]. The validity of this result is questionable 
based on the reasoning given below [322]. 

In the model of Refs. [300,321], the width of the front W increases with 
the external load approximately as a power law, and eventually saturates. As 
in gradient percolation, the saturated width W* scales with the gradient of 
the damage profile (see Eq. (jHZI))- Since in Ref. [300] ^ ~ L, where is the 
lattice size parallel to the front, one can combine the initial dynamic scaling 
with the saturated width into a "Family- Vicsek" -like scaling form W{Lx,t) = 
L°'f{t/L^), and conclude that the fronts are self-affine interfaces with = a. 
In gradient percolation, however, a can not be interpreted as a roughness 
exponent, since as we discussed above the front is self-similar [322]. If we 
remove the overhangs from the front through a solid-on-solid projection, we 
obtain a single valued interface h(x) that can be analyzed using conventional 
methods. The result, shown in Fig. 12^1 is an artificial multiscaling for the 
correlation functions 

(|/i(x + x') -/i(x')r)^/'^ (88) 

where = 1 for q < 1 and Q = 1/q for ^ > 1 The multiscaling behavior is 
easily understandable: for q < 1 the scaling is sensitive to the fractal, hence 
= 1 as expected for isotropic scaling. For q > 1 the scaling is dominated by 
the big jumps, which are also apparent in Fig. |2H1 [323]. 

In summary, if the front is due to a gradient percolation mechanism then its 
geometry is not self-affine, even when the gradient scales with L as in some 
fracture models [321]. The same reasoning applies to the theoretical arguments 
put forward in Ref. [324] to explain the roughness of the cracks in the random 
fuse models. The scaling theory of Ref. [324] rests on a number of assumptions: 

(i) in the limit of strong disorder fracture occurs at a percolation (possibly 
correlated) critical point; 

(ii) for weak disorder crack localization occurs, producing a quadratic damage 
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Figure 38. A fractal front in gradient percolation. We also display the solid-on-solid projection of 
the front which is single valued but not self-affine. Courtesy of A. Baldassarri, see also [323]. 

profile; 

(iii) the width of the damage profile is a linear function of the lattice size; and 

(iv) weak and strong disorder are characterized by the same v exponent. 

Given these assumptions, the authors of Ref. [324] argue that in the localiza- 
tion regime the crack is formed as a front in gradient percolation, leading to a 
scaling relation a = 2i//(l + 2u), where the factor 2 (not present in Eq. (|87j) ) 
comes from assumptions 2 and 3. As discussed above, the identification of a 
with is incorrect. In addition, the analysis of large scale numerical simu- 
lations yields a series of additional problems: (i) the existence of a reliable 
percolation scaling is dubious; (ii) the damage profile is not quadratic; (iii) 
the width of the damage profile is not a linear function of L [325] . 



6 Numerical simulations 

Here we review the main results, obtained from numerical simulations of sta- 
tistical models for fracture for the kind of models discussed in Section [l] We 
present a detailed discussion of the algorithms employed for such simulations 
in the appendix. The main work horse for these studies has traditionally been 
the RFM, which therefore will garner the most attention in what follows. 
The main motif of this section is to summarize the phenomenology of the re- 
sults, from the RFM and similar models, and to emphasize the relationship 
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Figure 39. The high order correlation function for the soUd-on-solid projection of a gradient 
percolation front. Since the original front is self-similar, the projection displays an artificial 
multi-affinity with effective exponents = 1/g for q > 1. Courtesy of A. Baldassarri, see also [323]. 

with the fracture mechanics framework outhned in Sec. |21 and the theoretical 
approaches reported in Sec. [3 For simphcity we defer a general comparison 
between simulations and experiments discussed in Sec. |31to the conclusions. 

We first discuss the scaling of the stress-strain characteristics (Sec. 16. If) , 
and its interpretation in terms of damage mechanics (Sec. One issue is 
whether the numerical data can be interpreted in terms of a scaling picture or 
framework from statistical physics. Damage as a thermodynamic quantity is 
then further explored analyzing the spatial distributions of local damage, the 
development of localization and internal correlations. The two central ques- 
tions here concern the nature of the approach to peak load, what damage can 
tell about it, and the role of disorder in the process. 

Next, we discuss the topics that are more relevant for statistical experiments: 
fracture strength, crack morphology and avalanche precursors. In particular 

(i) We analyze the properties of the strength distribution and the associated size 
effects, including a discussion of the effect of disorder in notched specimens. 
Such results are of importance due to the relation to weakest-link theory 
and due to the engineering applications fSec. 16.3)) . 

(ii) The surface properties of cracks are considered from the viewpoint of estab- 
lishing the true scaling indicated by the data: what is the minimal set of 
exponents that suffices, and what is the role of disorder and dimensionality 
rSeclOl). 

(iii) The avalanches in fracture models are discussed based on general scaling re- 
lations for systems with crackling noise. Again, it is important to understand 
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how the specific scenario (i.e., type of interactions, disorder, dimensionahty) 
influence the avalanche behavior fSec. I6.5)l . 

6.1 The I-V characteristics and the damage variable 

The first basic information that can be obtained from numerical simulation 
is the constitutive response of the model, in the form of a stress-strain curve 
for tensorial models or, more simply, of a current- voltage (I — V) curve for 
the RFM. This quantity is particularly important in connection with damage 
mechanics (see Sec. 12. 5|) . since lattice models allow to test possible relations 
between the microscopic damage state and the macroscopic constitutive be- 
havior. Here for the sake of clarity we focus on the scalar case (i.e. the RFM) 
where the damage variable D is simply defined. The same discussion can be 
translated to generic (tensorial) fracture models considering the stress strain 
characteristics. These models have in general the important feature, that they 
allow in contrast to mean-field -like approaches the existence of damage cor- 
relations, which can thus be explored. 

Figure EHl presents a typical I — V response, from a 2D RFM (triangular 
lattice), with a system size L = 64. One can see damage accumulation, and 
then signatures of crack growth which becomes unstable after the peak current. 
While initially the stress concentrations are masked by the disorder, in the last 
stage one of the microcracks wins over the others and dominates. 

To obtain a faithful representation of the macroscopic response one needs 
to perform an ensemble average over many sample configurations. This can 
naturally be done fixing three different control parameters: current, voltage, or 
accumulated damage (i.e. number of broken bonds). Appropriate rescaling of 
the curves is then necessary to obtain intensive quantities (i.e. independent of 
the system size). For this reason the voltage V is the most natural parameter, 
since the current excludes the post-peak load part of the I — V response, and 
the damage is - as is discussed below - not an intensive variable, at least in 
the post peak regime. 

An old idea, inspired from self-affine interface growth, was to attempt a 
Family- Vicsek (FV) [326] scaling for the I — V curve (see also Ref . [327] ) . This 
entails the scaling form 



The scaling function (/>(x) has to fulfill the condition that for small V the 
response is independent of system size. There is an important detail, not fully 
appreciated in the literature: usually in the context of surface growth the use 
of Eq. (|5^ implies that cj) const > for x 3> 1. In fracture, however. 
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Figure 40. Typical I — V response of a 2D RFM triangular lattice system of size L = 64. The thin 
line is the evolution that follows a quasistatic dynamics in which bonds are broken one at a time. 
The envelope of the I — V response (thick line) represents the voltage controlled response. 

the current becomes zero for large enough V: in that Hmit (j) vanishes and is 
thus not monotonic. For small x the function has to scale as power-law of its 
argument, which directly implies that the FV collapse is not applicable unless 
the peak load has a power-law scaling form {Imax ~ , for some y), and the 
same holds also for the corresponding voltage. In addition, the scaling of (f) 
couples the exponents and Q.2'- since for small V , I (x V , it follows that 
(p ^ X and = 0.2. While this discussion is valid in general, for arbitrary 
models and experiments, we consider as an example the numerical data for 
the RFM. 

Figure lUT a) presents the scaling oi I — V response based on the FV scaling. 
The plots in Fig. HlT a) correspond to Q.i = 0.2 = 0.75 as proposed in Ref. [326]. 
On the other hand, an alternative ad-hoc scaling form 



{logLY 



(90) 



has been proposed by Sahimi and collaborators, in Refs. [328,329]. The values 
indicated are ~ litO.l and V' — 0.1. Figure ETT b) presents the corresponding 
scaling oi I—V curves based on Eq. (I^U)) . This suggests that scaling form based 
on Eq. H90() is de facto more appropriate than the FV one, but of course one 
notes that the scaling is far from perfect. The logarithmic factor in Eq. H90() 
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arises from an attempt to correctly take into account the strength size effect. 
For such a scahng to hold the scaling function (j)2 would also need to contain 
logarithmic terms so that the small-voltage limit is independent of L. Their 
origin would be unknown. Notice that the analysis presented here is also valid 
for experiments, and for instance questions the application of FV scaling to 
the experiment of Otomar et al. [177]. 




Notice that for the data in Figs. |^ ensemble averaging is performed us- 
ing the accumulated damage as a control parameter as originally done in 
Refs. [326,328-330]. This produces an unphysical hysteresis in the curves, 
which corresponds to unstable crack growth in voltage or current controlled 
cases. This effect was not apparent for the small lattices sizes at hand in 
Refs. [326,328-330], but becomes clear for larger lattices employed here. To 
avoid this problem, it is more natural to average the envelopes of the I — V 
curves, using the voltage as a control parameter. The total current I flowing 
through the lattice in the RFM is given hy I = T, V where S is the total 
conductance of the lattice system expressed as 



(91) 



Y2 

and aij is the local conductivity of the fuse between nodes i and j. A simple 
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normalized form of the I — V curve may then be obtained as 



IL 



V 

(92) 



Nel Nel L' 

where Sij denotes the enhancement of current (stress concentration factor) due 
to fuse burnouts and ^ is voltage drop per unit length. It should be noted that 
Sij = for all the dangling or dead ends (using percolation theory language) 
that do not carry any current (stress) even though they are not broken. 

Figure presents the ensemble averaged I — V response normalized as in 
Eq. ()92() . showing a perfect collapse up to peak load, with strong deviations 
thereafter. One can easily extract from this curve the behavior of the damage 
variable D(y) = 1 — S/Sq, where Eg is the conductivity of the intact lattice. 
Figure 1331 shows that damage variable D(y) grows linearly and remains an 
intensive variable up to the peak load, beyond which size dependent effects 
become apparent. Similar damage variable definitions are often used in the 
engineering literature [25]. It is clear that they are suitable to represent the 
extent of damage only within the early stages of damage evolution, where 
the damage levels are moderate and the damage is not localized, but become 
clearly inadequate around the peak load and especially in the softening regime 
of the stress-strain curve. In these regimes, D is strongly size-dependent and 
hence cannot be considered as a thermodynamics based intensive damage state 
variable. Of course, one should keep in mind that due to the "size-effect" 
asymptotically Imax{L) may go to zero anyway, which implies that Diy) would 
not exist in a well-defined sense for any particular value of V . 

In order to better understand the behavior of the damage variable, it is in- 
structive to observe the degradation of the conductivity of a single sample as a 
function of the fraction of broken bonds pb- As is apparent looking at Fig. I44[ 
the conductivity before the peak load follows nicely the prediction of effective- 
medium theory, considering the contribution to ph coming both from broken 
and dangling bonds. This means that there are no strong sample-to-sample 
variations (nor as we discuss below a strong L-dependence). The agreement 
with the conjecture implies that damage is spread on large scales homoge- 
neously in the sample with no localization or strong correlations. Hence the 
definition of a RVE is possible and a continuum damage description is appro- 
priate. The lack of a "critical scenario" with a divergent correlation length ^ 
in any case implies that the scalar damage variable D could be simply writ- 
ten as a power series in the control parameter 1/, D = d^V + di^^ . . ., with 
I = (1 — D{y))TiQV . Then, the study of damage becomes that of the coeffi- 
cients di and their dependence on model, dimensionality, and disorder. 

In the post peak regime, however, the conductivity deviates abruptly from 



96 



the effective medium prediction, with strong sample-to-sample and lattice size 
dependent fluctuations. This brings back the question of ensemble averaging: 
by averaging together several curves as the one reported in Fig. 1441 the abrupt 
character of the discontinuity is partly lost and one obtains smooth curves as 
the ones reported in Fig. I45r a). This effect becomes less and less important 
as the sample size is increased, but has lead to some confusion in the early 
literature, when only small samples were available. Alternatively, the averaging 
of the conductivity degradation responses may be performed by first shifting 
the data by the respective peak load locations {pp,T,{pp)). In this way, one 
preserves the main character of the curves as shown in Fig. I45r b). 
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Figure 42. Normalized ensemble averaged I — V response based on Eq. 1921 for different 2D RFM 
triangular lattice systems of sizes L = {4, 8, 16, 24, 32, 64, 128, 256, 512}. 



6.2 Damage Distribution 

Fuse networks can be used to study the microscopic dynamics of damage. In the 
Sec. 16.11 we discussed the macroscopic characteristics of the RFM in terms of 
the I-V curves. The increase in the damage variable D (in I = {l—D(y))'SoV)) 
is concomitant to the growth of correlations due to microcracking controlled 
by the distribution of quenched disorder. When the disorder is narrowly dis- 
tributed (weak disorder), it is known that damage is localized and the network 
breaks down without significant precursors [165]. For infinitely strong disorder, 
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Figure 43. Ensemble averaged damage variable D(V) for different 2D RFM triangular lattice 
systems of sizes L = {4, 8, 16, 24, 32, 64, 128, 256, 512}. After peak load the damage variable is not 
intensive. Also, the nonlinearity in the plots of D{V) close to 1.0 is an artifact of I — V averaging 

close to fracture. 
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Figure 44. Typical conductivity (stiffness) degradation in various 2D RFM samples: triangular 
lattices of sizes L = 256 (dotted), 512 (solid), 1024 (thin dashed). The thick dashed line 
corresponds to the degradation of the conductivity predicted by effective medium theory. 
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Figure 45. Average conductivity scaling in a 2D RFM triangular lattice system, (a) Averaging 
performed at constant raj, number of broken bonds, (b) Averaging by first shifting each of the 
sample responses by their respective peak load positions. 

the damage accumulation process can exactly be mapped onto an uncorrelated 
percolation problem [314]. The interesting situations are to be found for strong 
but finite disorder, where a substantial amount of damage is accumulated prior 
to failure. In these conditions, one can try to study several quantities such as 
the scaling of damage at maximum current, or the geometric characteristics of 
the ensembles of broken fuses along I — V curves. One interesting question is 
in addition related to the fluctuations of damage at the I — V curve maximum. 

Figure |3ni presents the snapshots of damage evolution in a typical 2D RFM 
simulation of size L = 512 and uniform disorder distribution in [0, 1]. From this 
figure, at least visually, it is apparent that damage evolves in a fashion that 
does not foretell the final crack up to the peak load (see Fig. 146b ). beyond which 
the localization of damage occurs. Next we discuss various results concerning 
the distributions and fluctuations of damage, both at peak load and at final 
failure. 



6.2.1 Scaling of Damage Density. We first consider the evolution of the 
number of broken bonds at peak load, Up, which excludes the bonds broken 
in the last catastrophic event, and the number of broken bonds at failure, 
Uf. The quantity Up is of interest since it tells how the damage is distributed 
before failure occurs. If the damage were to be localized along a single row, 
the number of broken bonds at peak load and failure would be 0(1) and 0{L), 
respectively, whereas screened percolation would result in nj ~ N(,i, where Ni,i 
is the total number of bonds in the intact lattice. 
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Figure 46. Snapshots of damage in a typical 2D triangular lattice system of size L = 512. Number 
of broken bonds at the peak load and at failure are 83995 and 89100, respectively, (a)-(i) represent 
the snapshots of damage after breaking ni, number of bonds, (a) nt = 25000 (b) ni, = 50000 (c) 
rib = 75000 (d) rit = 80000 (e) = 83995 (peak load) (f) rit = 86000 (g) = 87000 (h) 
rib = 88000 (i) nt = 89100 (failure) 



Figure E7I presents Up as a function of the lattice size N^^i for 2D RFM and 
for 2D RSM with uniform threshold disorder. The data displays a reasonable 
power law behavior Hp ~ N^i, with b = 0.91 - 0.93 [182,325,326,331,332]. In 
other words, the damage is sub-extensive reflecting a size-effect at the peak 
load that decreases with N^i. The difference between the RFM and RSM mod- 
els is marginal. Some systematic deviations from the scaling form Up ~ Nj^i 
can perhaps be seen by plotting Up/Nj^i vs N^i (see Ref. [325]). Such attempts 
would indicate that the damage data be equally well fit by a linear law times 
a logarithmic correction Up ~ N^i/ \og{Nf,i) as suggested in Ref. [333]. This 
also implies that in the thermodynamic limit the peak load damage vanishes. 
For strong disorder (power law distributed with A = 20, see Sec. 14. 1() the data 
depicted in Fig. 1471 implies a volume-like 6 ~ 1, indicating behavior similar to 
the percolation limit in the range of L used. The three dimensional results for 
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the RFM are similar [334,335], with Up ~ N^^ and b = 0.72 for weak disorder 
(a uniform distribution in [1/2,3/2] [335]) and b = 0.99 for stronger disorder 
(a uniform distribution in [0, 1] [334]). The weak disorder scahng can also be 
recast as 2± corresponding to a localization along a plane. 




Figure 47. Scaling of number of broken bonds at peak load for two-dimensional triangular and 
diamond random fuse lattices, and triangular spring lattices. In all these models disorder is 
uniformly distributed in [0, 1] and the scaling exponents are very close to each other. The last curve 
represents data obtained using a power low disorder distribution: the exponent changes to 6 = 1. 



The other interesting quantity is the average number of broken bonds in the 
post-peak regime, (n/ — Up). This corresponds to the size of the last catas- 
trophic avalanche, which is trivially so in particular in the current-driven case. 
It also exhibits a power law scaling N^^ with ^ = 0.72 in 2D RFM triangular 
lattices and 9 = 0.69 in 2D RFM diamond lattices [325], while 9 = 0.68 in 2D 
RSM triangular lattices [182], and 9 = 0.81 in 3D RFM cubic lattices [334]. 
If instead one uses again the A = 20 case depicted above, the scaling with 
Nf,i seems to be linear (9 = 1), i.e., extensive similar to Hp [325] whereas for 
weaker disorder one has 9 < b. 

If n J scales extensively, one may then ask what is the finite-size correction to 
an asymptotic density defined as pc = limi_^oo nf/N^i, in analogy to percola- 
tion as argued by Hansen and co-workers [324,336]. A percolation-like scaling 
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for the finite-L fraction of broken bonds would be given by 

Pf-Pc^L-T. ^N',^" (93) 

Here, Pf = and z/ would define an exponent similar to that of the correla- 
tion length exponent in percolation. While small system size simulations could 
be consistent with Eq. (pS)) . Ref. [325] shows that percolation scaling can be 
excluded and there is no hidden divergent lengthscale defining an exponent 
v. This is in agreement with the first-order, nucleation picture of fracture 
that one gains from mean- field -like models (Sec. 15. ill . A similar conclusion is 
corroborated by considering other measures related to a characteristic length- 
scale (which may be called a correlation length or a "RVE size" , as discussed 
in Sec. 13) [331,333]. 

To study further the question of damage fluctuations at peak load and fail- 
ure we next consider the distributions ensemble averaged over samples. The 
cumulative probability distribution for the damage density at the peak load is 
defined as the probability Ilp{pb, L) that a system of size L reaches peak load 
when the fraction of broken bonds equals ph = where Ub is the number 
of broken bonds. Scaling such distributions using the average and standard 
deviation one obtains a good collapse of the distributions (Fig. Ii8|) for both 
random fuse and random spring models in d = 2 and it seems that this behav- 
ior is independent of lattice structure and disorder [325]. Similar results are 
found in three dimensions. 

The cumulative probability distribution for the damage density at failure is 
defined similarly, nj(pb, L) [325], and there seems to be no essential difference 
between the scaled forms of IIj and lip. This is in contrast with some claims 
in the literature wherein a percolation scaling nj(p6,L) = H{(pb — Pc)L^^'^) 
was suggested for Yif{pb,L), and the width of the distribution is assumed to 
originate from a correlation length-like scaling [324] . In the random fuse model, 
the scaling forms of cumulative distributions of damage density at peak load 
and failure seem to follow the same scaling function, Il{p) = Ilf{pf) = Ilp{pp) 
[325]. This does not seem to hold for spring lattice systems [182], possibly due 
to the non-local character of the loss of rigidity in central-force models (Sec. 
14. 2|) . We note finally, that the scaling functions display Gaussian forms in the 
case of the RFM (see the inset of Fig. 0H1 and Ref. [182] for the peak load 
of 2D RSM). The Gaussian distribution is consistent with the absence of a 
divergent correlation length at failure [331,333], which can be understood by 
observing that the damage density pj is the sum of local random variables (the 
local damage, Sj), each originating from the contribution of a single RVE. The 
central limit theorem yields then a Gaussian distribution, only if correlations 
in the field Sj are short-ranged as is the case if a RVE can be defined. 
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Figure 48. Comparison between the cumulative probability distributions of the fraction of broken 
bonds at the peak load is presented for the RSM and RFM. For the RSM, triangular lattices of 

sizes (L = 8, 16, 24, 32, 64, 128, 256), and for the RFM, triangular lattices of sizes 
(L = 16, 24, 32, 64, 128, 256, 512) are plotted. The normal distribution fit in the re-parameterized 
form is presented in the inset. In the RFM case, collapse of cumulative probability distributions at 
the peak load for different lattice topologies (triangular and diamond) and different disorder 
distributions (uniform and power law) is presented in Ref. [325]. 

6.2.2 Damage Localization. The question of how damage locahzes in lat- 
tice models is a problem with many facets. An important question is how 
the damage variable D is connected to Uh/N^i and to the effective RVE size. 
One has to separate here two issues, the first of which is the development of 
a macroscopic damage profile in a sample, and the second being how local 
microfailures become correlated. The first issue can be resolved by evaluating 
the disorder averaged damage profiles from simulations. However, a posteriori 
analysis of damage profiles is tricky since the final crack position fiuctuates, 
i.e, it can occur anywhere along the loading direction of the lattice system, 
and the fact that the final crack position is in general surrounded by a damage 
"cloud" . To overcome this problem, one can use different techniques to deter- 
mine the location of the "crack" along the macroscopic current fiow direction. 
Some of these techniques are based on shifting individual damage profiles by 
the center of mass of the damage, or averaging the power spectrum (square of 
magnitude of the Fourier transforms) of individual damage profiles [325]. 

Figure 011 presents the average damage profiles for the damage accumulated 
up to the peak load by first shifting the damage profiles by the center of 
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mass of the damage and then averaging over different samples. This result, 
displayed in Fig. 0^1 for 2D RFM model with uniform disorder distribution, 
clearly shows that there is no localization at peak load. Similar behavior has 
been observed in various other cases; namely, the 2D RFM simulations with 
power law threshold distribution (A = 20) [325] , the 2D RSM simulations with 
uniform thresholds distribution [182], and 3D RFM simulations with uniform 
thresholds distribution. For small system sizes, L, there are residual effects 
due to the non-periodicity of the lattice system in the loading direction. All 
in all, such profiles suffice to show that the localization profiles measured at 
failure in the literature [324,333,336,337] are due to the damage accumulated 
between the peak load and failure. Thus, in the following, we consider the 
damage profiles Ap(y) for the damage accumulated between the peak load 
and failure. Figure EOl presents a data collapse of the average profiles, using 
the form 

{Ap{y, L))/(Ap(0)) = f{\y - L/2\/0, (94) 

where the damage peak scales as (Ap(0)) = L~^, with x = 0-3 and the 
profile width scales as ^ ~ L", with a = 0.8 (see Fig. I50|). Similar results are 
obtained for the 2D RSM [182] and 3D RFM [338] models: the profiles display 
exponential tails in both cases, with the exponents for the peak being x = 0.37 
and X = 0-15 and for the profile width a = 0.65 and a = 0.5, respectively. 
The Q-exponents clearly have to do with the final crack roughness, which is 
discussed in Sec. 16.41 It is an open question as to how the dynamics in this 
catastrophic phase should be characterized, in order to explain the value of a 
(or possibly that of conversely), and more work is needed in this respect. 
This would also call for studies with simplified models (Sec. 15.1]) . 

One may also consider as a measure the "width of the damage cloud" , which 
is defined asW = {{{yb — yb)'^))^^'^- Here yb is the y coordinate of a broken bond 
and a disorder average is performed. At peak load, W ~ L/y^l2) for both 
2D RFM and RSM lattice systems [182,325,338], consistent with uniformly 
distributed damage profiles at the peak load. In the post-peak regime, if the 
average is restricted only to the bonds broken in the last catastrophic failure, 
one has W ~ L^-^^ for 2D RFM [325], and W ~ L^-^ for the 3D RFM. These 
non-trivial scaling exponents in the post-peak regime indicate that the damage 
profiles exhibit a localized behavior in the post-peak regime in agreement with 
the results presented above. 

The presence of an abrupt localization can also be seen by considering the 
difference in the accumulated damage between the window containing the final 
crack (and damage) and any other arbitrary window. Let ni and denote 
the associated numbers of broken bonds. Figure |^ presents the evolution 
of (n; — Ua) with damage accumulation for 2D RFM lattice system of size 
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Figure 49. Average damage profiles at peak load obtained by first centering the data around the 
center of mass of the damage and then averaging over different samples. Data are for RFM 
simulations with uniform thresholds disorder. 

L = 512, averaged over 200 samples. This result shows that (n; — Ua) ~ up 
to the peak load, and a steep increase in {ni — Ua) in the post-peak regime 
indicating abrupt localization. 



6.2.3 Crack clusters and damage correlations. An analysis of the clus- 
ters of connected broken bonds is the easiest way to consider the correlation 
properties of damage. This topic is also of interest due to the relation between 
the largest defect cluster and the fracture strength distributions discussed in 
Sec. 12.41 [9,329,339,340]. In the RFM with strong disorder damage accumu- 
lation is non-trivial. One could postulate that the size-scaling and strength 
distribution correspond to the statistics of the largest crack present in the sys- 
tem at peak load, while crack represents the result of damage accumulation 
prior to the peak load. 

The relevant questions are: (i) what kind of distribution is followed by the 
crack clusters in the strong disorder case just before the appearance of an 
unstable crack, and (ii) whether there exists a relationship between the dis- 
tribution of maximum cluster size at peak load and the fracture strength 
distribution? 

Figure presents the connected cluster size distribution at the peak load 
in a 3D RFM for various L. The plots indicate that simple power-law and 
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Figure 50. Data collapse of the localization profiles for the damage accumulated between peak 
load and failure. The profiles show exponential tails. Data are for RFM simulations with uniform 

thresholds disorder. 
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Figure 51. Localization of damage (rt; — rta) in a 2D RFM simulation of lattice system size L = 512 
averaged over 200 samples. Localization of damage occurs abruptly at the peak load (ni,/np = 1). 
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exponential representations are not adequate, and a similar result is obtained 
in 2D as well [341]. For completeness, one should note that the distribution 
of the largest crack cluster at the peak load seems to follow a Lognormal dis- 
tribution as demonstrated in Fig. |S31for various system sizes. It is interesting 
to note that the fracture strength data presented in the following subsection 
indicate that a Lognormal distribution represents an adequate fit for the frac- 
ture strength distribution in broadly disordered materials. At first glance, this 
result appears to be trivial: assuming that the stress concentration around the 
largest crack cluster is given by a power law, the fracture strength follows a 
Lognormal distribution since the power of a Lognormal distribution is also a 
Lognormal distribution. However, as noted in Ref. [341], a close examination 
of the largest crack cluster size data at peak load and the fracture strength 
data reveals that the largest crack cluster size and the fracture strength are 
uncorrelated [341]. Indeed, quite often, the final spanning crack is formed not 
due to the propagation of the largest crack at the peak load, but instead due 
to coalescence of smaller cracks (see Figs. 6 and 7 of Ref. [341]). 




In(s) 

Figure 52. Damage cluster distribution at the peak load in cubic lattices of system sizes L = {10, 
16, 24, 32, 48} in log-log and log-linear scales (inset). The cluster distribution data for different 
lattice system sizes is neither a power-law or an exponential distribution. 

As a further method to assess the presence of correlations, we consider the 
ensemble averaged response of entropy, which has also been used in experi- 
ments (see Fig. 1^ in Sec.|3|. The entropy is defined as 



S = -(^Qi log{qi)), 

i 



(95) 
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Figure 53. Re-parameterized form of Lognormal fit for the cumulative probability distribution of 
largest crack cluster size at the peak load in cubic lattices of system sizes L = {10, 16, 24, 32, 48}. 
fip^ and refer to the mean and the standard deviation of the logarithm of largest cluster 

sizes, Sp, at peak load. 



where qi is the fraction of broken bonds within the i^^ window (or section). 
The lattice system is divided into i = 1, . . . ,mf,o^ independent measurement 
windows. 

A useful choice is to take the windows to be of width Ay, along the current 
flow direction [333] . We now present data for the load-displacement response of 
a typical RFM simulation, divided into 12 segments, with six equal segments 
each before and after the peak load. Let Dj = 1,2, . . . , 12 refer to each of 
these segments. The entropy S of the damage cloud is calculated at the end 
of each of these Dj stages using different section widths. For each system 
size L, we considered A^^^ox number of different section widths, where Nh^r^ = 
^"logi"^^ ^' the different window sizes is denoted by a box size index 

i?/ = (0, 1, • • • , Nijox — 1)) with the corresponding section widths Bsize = 2^^. 
Thus, there are mi,ox = number of sections (boxes) over which is 
performed (i.e., i = 1, . . . , nibox)- 

Fig. presents the entropy S versus Dj, calculated using different section 
widths for a uniform threshold distribution (A = 1). For damage accumu- 
lation up to the peak load, the entropy increases and attains a maximum 
{Smax ~ log{mi,ox)) corresponding to a completely random spatial damage 
distribution, followed at the peak load by a localization phase in which the 
entropy decreases. In the inset of Fig. [S^the entropy at the peak load versus 
the box size index Bj is plotted for uniform and power law threshold distribu- 
tions, for system sizes L = (64, 128, 256, 512, 1024). The slope of the graph in 
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the inset equals —log{2) = —0.693, which matches exactly with the slope of the 
maximum entropy Smax versus Bj since Smax = log{nn,ox) = log{L/Bsize) = 
log{L/2B') = log{L)-Bilog{2). 

1 .0005 1 ~ ~ , ~ ~ 1 




2 4 6 8 10 12 

Damaged Stage 

Figure 54. Entropy S versus the damaged stage index Dj for a triangular lattice of size L = 512 
with uniform threshold distribution (A = 1) threshold distribution. Symbols are 
Bi = (0, 1, ■ ■ ■ , Niax ~ 1) = {*! +1 °i The inset presents the entropy at the peak load 

versus the box size index for uniform and power law threshold distributions for 
L = (64, 128, 256, 512, 1024). The data for smaller systems is shifted such that the points 
correspond to the same relative box size. The slope of the graph in the inset is equal to 
—log{2) = —0.693 corresponding to randomly distributed damage. 



6.3 Fracture strength 

6.3.1 The fracture strength distribution. Before discussing the results of 
numerical simulations let us recall the main theoretical points of the strength 
of disordered materials. As discussed in Sec. 12.41 the widely used Weibull 
and modified Gumbel distributions naturally arise from the extreme-value 
statistics if one assumes a particular form of the defect cluster size distribution 
in a randomly diluted network: (1) defect clusters are independent of each 
other, i.e., they do not interact with one another; (2) system failure is governed 
by the "weakest-link" hypothesis, and (3) there exists a critical defect cluster 
size below which the system does not fail, and it is possible to relate the 
critical size of a defect cluster to the material strength. In particular, if the 
defect cluster size distribution is described by a power-law, then the fracture 
strength obeys Weibull distribution, whereas an exponential defect cluster 
size distribution leads to the Gumbel distribution for fracture strengths [9, 
329,339,340]. Theoretical arguments along these lines have been put forward 
in the context of lattice models, as discussed in Sec. 15.41 A different theoretical 
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Figure 55. Examples of how the Gumbel and WeibuU forms fit into the percolation disorder RFM 

data (from Ref. [316]). 

approach is based fiber bundle models which yields in the case of ELS [14], a 
Gaussian strength distribution with variance proportional to 1/iV, where is 
the number of fibers [14]. The more interesting case of LLS [256], yields a mean 
fracture strength /// = (o"/) decreasing as l/log(A^), such that the strength 
of an infinite bundle system is zero. The fracture strength distribution follows 
Eq. (|59() , and is close to a Weibull form although it is difficult to determine its 
exact form [256]. Using renormalization techniques, an asymptotic size effect 
of the form fij ~ log(log(L)) was proposed in Refs. [339, 342-344] for LLS 
models. 

Duxbury et al. [22, 316, 317] studied the distribution of fracture thresholds 
in the randomly diluted fuse networks with bond percolation disorder (see 
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Sec. I5.4jl . Fig. 1551 shows simulation results [316] compared with the Gumbel- 
Ansatz. Subsequently, the fracture strength distribution was studied in a vari- 
ety of lattice models with different local behavior including central-force spring 
models [328,345], Born models [184], and bond-bending models [328]. These 
studies along with the analytical investigations based on largest stable defect 
size [165] concluded that the fracture threshold decreased with the system size 
L as a power of log(L) and the fracture threshold distribution is best described 
by a double exponential (modified Gumbel) distribution. In these randomly 
diluted disorder problems, the defect cluster size distribution is exponential 
far away from the percolation threshold and follows a power law close to the 
percolation threshold. Consequently, in general, a Gumbel distribution better 
fits the fracture strengths distribution far away from the percolation thresh- 
old and a Weibull distribution provides a better fit close to the percolation 
threshold [9,329,339,340]. 

In the case of RFM models with weak disorder, failure is dominated by the 
stress concentration effects around the crack tip, and hence the "weakest-link" 
hypothesis should be valid. Thus the fracture strength follows extreme-value 
theory and is distributed according to the Weibull or modified Gumbel distri- 
butions. The interesting scenario corresponds to the case of strong disorder, 
where neither Weibull nor Gumbel distribution applies well (see Ref. [346] for 
2D RFM, Ref. [182] for 2D RSM, and Ref. [334] for 3D RFM). There are two 
main reasons behind this behavior: 

(i) in the "weakest-link" hypothesis, the fracture strength is determined by 
the presence of few critical defect clusters, and is defined as the stress re- 
quired for breaking the very first bond in the system. In materials with 
broad disorder, the breaking of the very first bond ("weakest-link") does 
not usually lead to the entire system failure (as is manifested with the scal- 
ing Up ~ A'^^; with b = 0.93 for 2D RFM with strong disorder), and hence 
a fracture strength distribution based on the importance of very first bond 
failure may not be applicable. 

(ii) As discussed above in Sec. I6.2.3( the actual geometric crack cluster size 
distribution in the RFM at the peak load develops its own form, quite 
different from the initial defect cluster size distribution in randomly diluted 
networks. In general at the peak load of RFM simulations, the crack cluster 
distribution is represented neither by a power law nor by an exponential 
distribution (see Fig. |S21and Ref. [338]). Hence, neither the Weibull nor the 
Gumbel distributions fit the fracture strength distribution accurately. 

Figure EHl shows the fracture strength density distributions for random 
thresholds fuse and spring models using the standard Lognormal variable, 
^, defined as ^ = k!l!£lLJl^ where Uf refers to the fracture strength defined as 
the peak load divided by the system size L, and rj and C refer to the mean 



Ill 



and the standard deviation of the logarithm of cj/. The excellent collapse of 
the data for various fuse and spring lattices clearly indicates the universality 
of the fracture strength density distribution. 

In addition, the collapse of the data in Fig. ISBl suggests that P{a < Uf) = 
^{C)-, where P{a < aj) refers to the cumulative probability of fracture strength 
0" < (T/, ^' is a universal function such that < ^' < 1, and ^ = is 
the standard Lognormal variable. The inset of Fig. [SHI presents a Lognormal 
fit for fracture strengths, tested by plotting the inverse of the cumulative 
probability, <I>~^(P(cr/)), against the standard Lognormal variable, ^. In the 
above description, $(•) denotes the standard normal probability function. As 
discussed in Ref. [334], a Lognormal distribution is an adequate fit for 3D 
random fuse models as well. This further confirms the notion of universality 
of fracture strength distribution in broadly disordered materials. 
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Figure 56. Universality of fracture strength distribution in the random thresholds fuse and spring 
models. Data are for different lattice system sizes, L, corresponding to triangular fuse lattice with 
uniform disorder (L = {4, 8, 16, 24, 32, 64, 128}), diamond fuse lattice with uniform disorder 
(L = {4, 8, 16, 24, 32, 64, 128}), triangular fuse lattice with power law disorder 
(L = {4, 8, 16, 24, 32, 64, 128}), and triangular spring lattice with uniform disorder 
{L = {8, 16, 24, 32, 64, 128}). The inset reports a Lognormal fit for the corresponding cumulative 

distributions. 

The origin of a Lognormal strength distribution is not entirely clear from 
the theoretical point of view (for a heuristic argument see Ref. [346]). It is 
worth remembering that what determines the survival of a system is a quantity 
with strong sample to sample variations: the microcrack with the largest stress 
enhancement. In contrast to percolation disorder, this critical crack grows from 
the no-damage state in these simulations. In Sec. 16.2.31 we have discussed also 
the scaling of the connected crack clusters. It is worth remembering that the 
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stability of a given configuration results from two contradictory effects: stress 
enhancements and configurations (of fuses in the RFM) that are able to arrest 
the crack and make it stable (as is the case in LLS fiber bundles [276]). In 
experiments on paper, there are results that point out in a similar vein the 
combined role of local strength and local stresses, which can vary not only 
because of microcracking but also because of structural inhomogeneity [168]. 

6.3.2 Size effects. Two forms have been proposed in the literature to ex- 
press the size effects on mean fracture strength in randomly diluted lattice 
systems (percolation disorder): a simple power law scaling [9], consistent with 
the Weibull statistics, 

2 

/i/ ~ L~ ™ , (96) 

or a more complicated form [22,316,317,328,329,339,345,347], consistent with 
a modified Gumbel description, 

where rh is the Weibull exponent, and 6, A\ = log ^2 and B\ are constants. 

The inset of Fig. [^71 presents the scaling of mean fracture strength, ^u/, based 
on power law scaling /ij ~ L ™ . From the nonlinearity of the plot, it is clear 
that the mean fracture strength does not follow a simple power law scaling. 
On the other hand, the scaling form /x^ = + Bi in L ™ay represent the 
mean fracture strengths reasonably well even though the Gumbel distribu- 
tion is inadequate to represent the cumulative fracture strength distribution. 
However, the fit results in coefficients that are extremely sensitive; namely, 
S = -2.45, Ai ^ -18.3 {A2 « 10"^) and Bi « 20 (see for instance, Fig. 5b of 
Ref. [346] for 2D RFM data) consistent with the divergent behavior reported 
in Refs. [22,345] for randomly diluted systems. 

Alternatively, Ref. [346] suggests a scaling form for the peak load (peak 
current), Ipeak, given by Ipeak = CqL'^ + Ci, where Cq and Ci are constants. 
Correspondingly, the mean fracture strength defined as /x/ = ^rr, is given by 
fif = Co + Tpn-, where d = 2 in 2D and d = 3 in 3D. Using this scaling 

law, the two-dimensional RFM and RSM fracture strength data are plotted as 
shown in Fig. 1571 The results presented in Fig. [27| indicate that the exponent 
a = 0.97 for RSM and a = 0.96 for RFM using both triangular and diamond 
lattice topologies. 

The results for three-dimensional random fuse models are similar. A power 
law scaling form /// ~ consistent with a Weibull distribution does not 
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represent an adequate fit (see Ref. [334]). On the other hand, a scahng of 
the form = + Bi in L ^'^sults in a very high value of S = —7.5 for the 
exponent and the high values for the coefficients Ai = —156380 and Bi = 
85050 indicating the divergence behavior similar to that observed in 2D RFM 
and in randomly diluted systems [22,345]. Alternatively, using the scaling form 
Ipeak = CqL'^ + Ci , the mean fracture strength data for 3D RFM results in a 
value of a = 1.95 [334]. 

Since a very small negative exponent {a — d+\) is equivalent to a logarithmic 
correction, i.e., for {d — 1 — a) « 1, L'^~'^+'^ ~ (In L)^^, an alternative 
expression for the mean fracture strength may be obtained a.s ^if = ^i^^^, 
where and c are constants that are related to the constants Co and Ci. This 
suggests that the mean fracture strength of the lattice system decreases very 
slowly with increasing lattice system size, and scales as /i/ ~ JhTLyF-i with 
ip ~ 0.15, for very large lattice systems. The same asymptotic scaling relation 
(i.e., /Uj (;n ) conjectured in Ref. [328,329] for broadly disordered 

materials. 



6.3.3 Strength of notched specimens. Above, we have investigated the 
fracture of unnotched specimens, but from an engineering point of view, 
scaling of fracture strength in pre-notched specimens is of significant inter- 
est [52,73,74, 161,348]. The size effect in this case is often given by a scaling 
of the form /xj oc a~^/^, where a = qq + cj is the effective crack size, oq is 
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the initial crack (notch) size, and cj is the FPZ size surrounding the crack 
tips [52, 73, 74, 161, 348]. As we have noted in the previous section, a loga- 
rithmic size effect is observed in 2D and 3D simulations using the random 
thresholds models (RFM and RSM) with both uniform (A = 1) and power 
law (A = 20) thresholds disorder. The difference between a~^/^ size effect in 
the engineering literature [52,73,74,161,348] and the typical logarithmic size 
effect in the statistical physics literature is due to two aspects: i) initial relative 
crack size ao/I/, and ii) disorder. 

A simulation of an initially cracked (or notched) lattice system is reported 
in Fig. |SS1 for a 2D RFM. It illustrates the possible role of disorder in deter- 
mining the ultimate sample strength and its scaling with L and/or ao/L. The 
scenario corresponds to some initial distributed damage followed by the crack 
roughening, which at the trivial level indicates the effect of crack arrest due 
to strong bonds, followed finally by the point of catastrophic failure. 

The size effects obtained in such notched specimens using the 2D RFM with 
uniform thresholds disorder are presented in Fig. [2^1 showing the scaling of 
fracture strength with system size for a fixed initial crack size, uq/L. The 
results can be fitted by a power law of the type 

1^1 f oc (98) 

where the scaling exponent m is significantly influenced by ao/L: for small 
ao/L values, m is very small, and is equivalent to a logarithmic correction as 
in Sec. 16.3.^ while for large ao/L values, m approaches 1/2 as predicted by 
LEFM [52]. The reason for this behavior is that in the small uq/L regime, 
fracture is dominated by disorder, whereas in the large uq/L regime, fracture 
is controlled by the initial crack. 

A weak disorder leads to the standard m = 1/2 scaling, while for strong 
disorder we have a logarithmic size effect. One way to understand the data - 
qualitatively - is to take note of the fact that for sure a typical sample can 
accumulate damage upto a number of broken bonds Hp for which the notch- 
induced crack is still stable. E.g. for a crack of size "oq = 1" this implies, that 
Up can be interpreted roughly as a density of broken fuses, for which the inte- 
gral of the failure threshold distribution corresponds to the typical strength of 
the system if the bulk damage is excluded. This, in turn, is given by 7r/4 times 
the mean of the threshold distributions times a numerical 0(1) factor, coming 
from the number of fuses most affected by the current enhancement (at most 
two). The damage will then produce a number of microcracks competing with 
the flaw, which is of size oq- Increasing the control parameter uq/L implies, 
that the damage Up will decrease since the current enhancement increases. 
Eventually, the notch will win over the other flaws, but this is a slow process. 
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Figure 58. Evolution of damage in a notched sample. Although diffuse damage is present, a 
wandering fracture line starts from the notch breaking the entire lattice. 
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Figure 59. Size effects in notched specimens with uniform thresholds disorder. The effective 
scaling exponent of mean fracture strength depends on the ratio ag/L. 



6.4 Crack roughness 

Lattice models for fracture have been used in the past to analyze the rough- 
ness of the crack surfaces in various geometries. The RFM has been numer- 
ically simulated in two [335, 347, 349, 350] and three dimensions [351, 352] 
using various types of disorder. Other lattice models have also been simu- 
lated for this purpose, including the beam model [353], the RSM and the 
Born model [216,354]. Interfacial (planar) cracks have been studied with the 
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RFM [298], the beam model [299], and a long-range scalar model [321]. Due to 
numerical limitations, earlier numerical simulations mostly considered global 
measurements for the crack roughness, while more recent results explored other 
quantities and the possibility of anomalous scaling. 

Let us first discuss the widely studied RFM [335,347,349-352]. As discussed 
in Sec. 14. II the plastic version of the RFM yields on minimum energy surfaces, 
which in two dimensions corresponds to the directed polymer, characterized by 
C = 2/3. The results of numerical simulations of the brittle RFM yield values 
close to C — 0.7 for a variety of disorder distributions [349] suggesting that 
indeed the fracture surface would result from minimum energy considerations 
[335] . High quality data for the surface roughness in d = 2 have been recently 
published in Ref. [350] ruling out this interpretation. The only caveat would 
be if the effective disorder would become correlated due the microcracking, 
which can change the ^ [347]. 

Several methods have been devised to characterize the roughness of an in- 
terface and their reliability has been tested against synthetic data [291]. If the 
interface is self-affine all the methods should yield the same result in the limit 
of large samples. In Fig. l6Ub we report the local width for triangular random 
fuse lattices for different sizes L. The curves for different system sizes do not 
overlap even for / <C L which would be the scenario when anomalous scaling is 
present (Sec. l2.6TT|) . The global width scales with an exponent ( = 0.83 it 0.02 
(C = 0.80 lb 0.02, obtained for diamond lattices [350]). On the other hand the 
local width increases with a smaller exponent, that can be estimated for the 
larger system sizes as Cioc — 0.7. A more precise value of the exponents is 
obtained from the power spectrum, which is expected to yield more reliable 
estimates. The data reported in Fig. HlOb are collapsed using C — Cioc = 0.13 
(C ~ Cioc = 0.1 for diamond lattices). A fit of the power law decay of the 
spectrum yields instead Czoc = 0.74 {( = 0.7 for diamond lattice), implying 
C = 0.87 {C = 0.8 for diamond lattice). The results are close to the real space 
estimates and we can attribute the differences to the bias associated to the 
methods employed [291]. Although the value of — Qoc is small, it is signifi- 
cantly larger than zero so that we would conclude that anomalous scaling is 
present. This implies in analogy to other theoretical scenarios that there is 
an external transverse lengthscale, increasing as a function of L. The question 
now becomes, what about the origins of such a scale? 

While the local exponent is close to the directed polymer value C, = 2/3, the 
global value is much higher. It is also interesting to remark that for directed 
polymers anomalous scaling is not expected to be present, even with noise 
which has a more complex character [32]. The result obtained in other lattice 
models are partly similar, posing the question whether there is a broad uni- 
versality class for fracture surfaces of disordered lattices. In particular, some 
measurements reported in the literature are C = 0.7 for a RFM with dilution 
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disorder [347], (" = 0.87 for the beam model in the hmit of strong disorder [353], 
Cioc = 0-66 in the Born model for a variety of parameters [216]. 

An additional tool to characterize the roughness is provided by the dis- 
tribution p(W) of the crack global width. This distribution has been mea- 
sured for various interfaces in models and experiments and typically rescales 
as [347,355,356] 

piW)=p{W/{W))/{W), (99) 

where (W) ~ L'' is the average global width. As shown in Fig. |^ the dis- 
tributions can be collapsed well using Eq. (|99() for diamond and triangular 
lattices [350], providing additional evidence for universality. Moreover, the 
scaling function follows a Lognormal distribution. 

The results of simulations in d = 3 are more difficult to interpret coherently 
since the numerical data is limited to small system sizes and fewer statistical 
samples. The higher computational cost of 3D simulations compared to 2D 
simulations restricts the consideration of 3D simulations to relatively small 
system sizes and smaller statistical sampling. It was debated whether the frac- 
ture roughness of the RFM could be described by a minimum energy surface, 
corresponding to C — 0.41. Simulations of the model with weak disorder lead 
to a similar value, but the exponent increases for stronger disorder [335,352]. 
Ref. [351] proposes an extrapolation leading to C = 0.62, but the numerical es- 
timates fluctuate considerably, depending on the method employed and on the 
considered disorder distribution. Recently performed extensive simulations of 
3D RFM have obtained similar results. In particular, an analysis of the width 
and the power spectrum data results in a value of C — 0.52, with possible 
anomalous scaling corresponding to a local value of Qoc — 0.42 (see Fig. I62|) . 
Finally, results of three dimensional simulations have also been reported for 
the Born model (on FCC lattices), yielding (loc — 0.5 from power spectrum 
and local width calculations. 

An interesting point to note is the consistency between the scaling of final 
crack width and the scaling of localization length ^ discussed in Sec. 16.2.21 
For the 2D RFM triangular lattice system, we have (" = 0.83 it 0.02, which is 
consistent with the scaling of the localization length ^ ~ L^'^ used in the data 
collapse of damage profiles in the post-peak regime. The center of mass shifting 
of the damage profiles for the purpose of averaging the damage profiles makes 
this consistency between the final crack width and the localization length ^ of 
the damage profiles self-evident. In addition, this agreement supports the idea 
that the post-peak damage profile is predominantly due to the localization 
produced by the catastrophic failure, which at the same time results in the 
formation of the final crack. Similar behavior is observed in three-dimensional 
RFM; namely, the scaling of final crack width ^ ~ 0.52 consistent with the 
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Figure 60. (a) The local width w{l) of the crack surface for different lattice sizes of 2D RFM 
triangular lattice systems plotted in log-log scale. A line with the local exponent Cioc = 0.7 is 
plotted for reference. The global width displays an exponent C > Cioc-(b) The corresponding power 
spectrum in log-log scale. The spectra for all of the different lattice sizes can be collapsed 

indicating anomalous scaling. 
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Figure 61. The distribution of the crack width is universal for diamond and triangular random 
fuse lattices. A fit with a log normal distribution is shown by a dashed line. 



scaling of localization length ^ ~ LP'^ [338] of the post-peak damage profiles. 
Interestingly, the same consistency in scaling between the final crack width 
and the localization length of the post-peak damage profiles is observed in the 
2D random spring models as well. The final crack width scales as ~ L*^'^^ 
(see Ref. [216]), and the localization length 5^ of the post-peak damage profiles 
scales as C ~ L^-^^ (see Fig. 6 of Ref. [182]). 



6.5 Avalanches 

The loading cm've of fracture models, such as the RFM, is characterized by 
avalanches of failure events. As discussed in Sec. l5.ll fiber bundle models repre- 
sent a useful idealization providing a qualitative description of the avalanche 
phenomenology and their statistics. On this basis one would expect to find 
a power law distribution with a cutoff changing with the applied load. An 
important theoretical question is whether the scaling exponents measured in 
lattice models are universal, i.e. they are the same in different models. On the 
experimental side, the main challenge is to relate the avalanche exponents to 
the acoustic emission statistics recorded in experiments. In both cases, one 
can ask what is the relation of the avalanches, and the underlying dynam- 
ics, to crack formation and damage mechanics. In lattice model simulations 
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Figure 62. The power spectrum of the crack surface of the three dimensional RFM. The spectrum 
exponent + 1 = 1.82, corresponds to (^jo^ = 0.41. The global width scales with = 0.52 (see 
inset), so that the power spectrum should be normalized by 2(^ — (^loc) — 0.2 

one typically defines the avalanche size s as in the FBM (see Sec- by 
counting the number of broken bonds that fail without any further increase 
in the applied load. Although this avalanche size is not directly equivalent to 
the acoustic emission energy, under some assumption, it is possible to relate 
the two distributions as discussed in the following. 

The avalanche statistics of the RFM has been reported for two dimensional 
[253,286,287,350] and three dimensional simulations [352,357]. While initially 
it was believed that FBM and RFM were in the same universality class [253, 
286, 287] with exponent r = 5/2, recent larger scale simulations displayed 
deviations from the expected behavior. The puzzling fact is, however, that the 
estimated value of the exponent r was found to be dependent on the lattice 
topology. This is contrary to any expectation, since in critical phenomena the 
exponents depend only on symmetries and conservation law, but not on the 
type of lattice or similar small scale details. An opposite situation is found in 
d = 3 where early simulations reported a strong deviation from mean-field, 
r ~ 2 < TpBM = 5/2 [352] while more recent results seem in agreement with 
this value [357]. Let us discuss all these results in detail. 

First, one can define avalanches increasing the voltage (deformation) or the 
current (force), but we do not expect strong differences for these two cases, 
so in the following we always refer to the latter. If one simply records the 
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distribution considering all the avalanches, one observes a power law decay 
culminating with a peak at large avalanche sizes. As in the FBM, the peak is 
due to the last catastrophic event which can thus be considered as an outlier 
and analyzed separately, avoiding possible bias in the exponent estimate (this 
problem is present in Ref. [358] reporting simulations of a central force model) . 
When the last avalanche is removed from the distribution the peak disappears 
and one is left in these models with a power law (see Ref. [350]). On the other 
hand, the size distribution of the last avalanche is Gaussian [182,350,359]. 

Next, as discussed in Sec. 15.1. J] one should distinguish between distributions 
sampled over the entire load history, and those recorded at a given load. In the 
first case, the integrated distribution has been described by a scaling form [350] 

P{s,L)=s~^g{s/L''), (100) 

where D represents the fractal dimension of the avalanches. To take into ac- 
count the different lattice geometries, it is convenient to express 2D scaling 
plots in terms of Nf.i rather than L 

P{s,N,i) = s-^g{s/N^/'). (101) 

A summary of the results for two-dimensional lattices are reported in Fig. 1631 
showing that in the RFM r > 5/2, with variations depending on the lattice 
topology, while the result reported for the RSM is instead very close to r = 
5/2. As for the avalanche fractal dimension, D varies between D = 1.10 and 
D = 1.18, indicating almost linear crack (which would correspond to D = 1). 
Note that the last avalanche has been studied earlier in the context of damage, 
and its geometrical properties in the context of crack roughness. The variation 
of D should be compared to that of C in the same systems. Three dimensional 
results are plotted in Fig. |M] for the RFM [357] . Again the agreement with 
FBM results is quite remarkable. 

The distribution of avalanche sizes sampled at different values of the current 
I is obtained by normalizing the current by its peak value Ic and dividing the 
I* = I/Ic axis into several bins. The distribution follows then a law of the 
type 

p{s,r) = s~^eM-s/s*),, (102) 
where the cutoff s* depends on /* and L, and can be conjectured to scale as 



(l-/*)i^L^ + C7' 



(103) 
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where C is a constant. Integrating Eq. (|1U2() we obtain 

P{s, L) ~ s-(^+") exp[-sC/L^]. (104) 

which imphes r = 7 + cr. We recall here that in the FBM 7 = 3/2 and a = 1 
(see Sec. I5.1.I|) . Fig. ESI reports a series of distributions obtained for the REM 
for different values of the currents [350]. The best fit for the scaling of the 
distribution moments as a function of I and L (Eq. (|1U3() ). and a direct fit on 
the curves suggest 7 ~ 1.9 — 2, a ~ 1.3 — 1.4 in d = 2 and 7 ~ 1.5, cr ~ 1 in 
d = 3. 

In summary, the avalanche distribution for REM is qualitatively close to 
the prediction of the EBM (i.e. mean- field theory). Apparently non- universal 
(lattice type dependent) corrections are present in d = 2, but the problem 
disappears in d = 3. This could be expected building on the analogy with 
critical phenomena, where exponents are expected to approach the mean-field 
values as the dimension is increased. Simulations of other models have been 
reported in the literature: the distribution obtained in Ref. [182] for the RSM 
is very close to the EBM result (see Eig. IH3j) . The avalanche distribution has 
been also reported for the Born model (r ~ 2 [360]), the bond-bending model 
(r ^ 2.5 [287]), the continuous damage model (r 1.3 [175]). 

To conclude, we discuss the issue of the relation between the AE distribution 
and and the avalanche size distribution. These two are, contrary to usual 
indications in the literature, not the same. The same holds even more clearly 
for the acoustic amplitude distribution. Acoustic waves were incorporated in 
a scalar model in Ref. [227], and the AE energy distribution was measured, 
with an exponent P ~ 1.7. This exponent can be related to r if we assume 
that the typical released AE energy of an avalanche of size s scales as e ~ s'', 
with 5 = 2 (i.e. the scaling relation is r = 1 -|- 6{P — 1)). This was verified 
in the model of Ref. [227], but also in the REM defining e as the dissipated 
electric energy during one event seems to agree [125]. 

7 Discussion and outlook 

We would like to conclude the present review by trying to answer two es- 
sential questions: (i) what can be learned from statistical models in order to 
reach a better understanding of fracture in quasi-brittle materials and inter- 
pret relevant experiments and (ii) what kind of developments can be foreseen 
to improve the current theoretical/modeling framework. The first question 
involves as well an assessment of the possible advantages of the statistical 
mechanics strategy with respect to traditional damage mechanics approaches. 
On the other hand, answering the second question will also help to clarify the 
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Figure 63. A summary of two-dimensional results for avalanche distributions integrated along the 
entire loading curve. The results are for the RFM on triangular and diamond lattices [350] and for 
the RSM [182]. Data for different lattice sizes have been collapsed according to Eq. IIOH . A line 
with exponent r = 5/2 is reported for reference. 



limitations of statistical lattice models and their role in the future. 

The macroscopic behavior of lattice models is well described using the ter- 
minology of phase transitions and the ensuing picture is conceptually rather 
simple: in the stress-controlled case crossing the peak-load implies a first-order 
transition to a "broken state" . This means that regardless of system size and 
disorder strength, the correlations in damage are local albeit with a growing 
characteristic length. This lengthscale does not diverge however as the peak- 
load is approached in contrast to that of a second-order transition. Such a 
picture has automatic consequences for the scaling of thermodynamic quan- 
tities as the conductivity or stiffness, damage and strength distributions. In 
particular, scaling laws inspired from critical phenomena typically fail to de- 
scribe the size dependence of these quantities. The statistical properties are 
instead well captured by relatively simple laws, such as Gaussian or Lognormal 
distributions, whose precise origin, however, is not always clear and remains a 
topic for future research. There are still interesting issues left, such as the most 
appropriate mathematical framework for measuring the damage anisotropy 
and correlations, or understanding the dynamics of microcrack populations. 
It would be important to carry out an extensive analysis for vectorial fracture 
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Figure 64. Data collapse of the integrated avalanche size distributions for the three dimensional 
RFM. The exponents used for the collapse are t = 2.5 and D = 1.5. A line with a slope r = 5/2 is 

reported in Ref. [350]. 

models, even though such quasistatic brittle fracture models are expected to 
have the same kind of "phase diagram" in the statistical mechanics language. 
Presumably allowing for several modes of element failure will also influence the 
right description of damage, summarized by the classical damage mechanics 
question whether it makes sense to use scalar damage parameter D{x). 

Next we would like to summarize with a few points the picture obtained from 
the theoretical approaches discussed in Sec. [3 and the numerical simulations 
reported in Sec. Inland compare the results with the experimental data analyzed 
in Sec. 01 An additional question that naturally follows from statistical models 
is the role of internal disorder for macroscopic damage laws used in engineering 
applications. 



7.1 Strength distribution and size effects 

One of the successes of statistical models of fracture has been to shed light 
into the scaling of strength with sample size and the associated distributions. 
There are indications that the asymptotics of extremal statistics are not the 
only alternative to explain the strength distribution, which results instead from 
cooperative phenomena or fluctuations in the damage accumulation. While for 
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Figure 65. The avalanche size distributions sampled over a small bin of the reduced current I* for 
a diamond lattice of size L = 128. The dashed line represents a fit according to Eq. 11021 with 

7 = 1.9 (from Ref. [350]). 

weak disorder the Weibull distribution provides an adequate fit to the data, for 
strong disorder we observe a Lognormal strength distribution. Whether this 
is a general result valid at all scales or the signature of a crossover behavior 
to Gumbel or Weibull scalings, still remains to be explored. 

A possible pathway to understand these results could stem from fiber bundle 
models (see Sec. 15. Under LLS conditions the model displays an extremal 
type distribution and a Logarithmic size effect. Conversely under ELS the 
distribution becomes Gaussian and size effects are absent. A clever way to in- 
terpolate between LLS and ELS was suggested by Hidalgo et al. [277] through a 
long-range rule in which the stress is transfered as 1/r"' . The crossover between 
LLS and ELS is identified with 7 ~ d, which corresponds with the exponent 
of the Green function obtained from lattice models, such as the RFM. Hence 
these models lie exactly at the boundary between LLS and ELS and the Log- 
normal distribution could then arise from the competition between local and 
global effects in a strongly disordered environment. Here, we would like to 
call for further basic science -oriented experiments. In the materials science 
-community it is customary to use the Weibull scaling to describe strength. 
However, it is clear that it is very rare to have a situation in which one would 
a priori expect that the conditions underlying the Weibull statistics are valid. 

Lattice models are also useful to analyze the effect on strength coming from 
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the interplay between elastic interactions and disorder in notched samples. 
For instance, the data displayed in Sec. l6.3.^ show the various mechanisms by 
which disorder can mask the classical LEFM prediction involving a I/^/oq- 
type scaling. The RFM has typically a rather small fracture process zone 
which decreases in size, for a fixed notch-to-sample size -ratio a^jL. The simple 
lesson we can draw from the model is that in practice one may, for a limited 
range of sample sizes, expect to get almost any strength exponent (in the 
law i7j oc ) between m = (or logarithmic behavior) and the m = 1/2 
prediction indicated by LEFM. 

It is interesting to close the present discussion by analyzing the predictions 
for strength coming from two other theoretical frameworks commonly used 
to interpret statistical properties of fracture, namely percolation (Sec. I5.4jl 
and depinning (Sec. 15.^1) . Percolation describes rigorously the infinite disorder 
limit of lattice models, where stress enhancement is irrelevant, the bonds fail 
one by one, and no real size effect is expected besides the small corrections 
associated with finite size scaling. The infinite disorder limit is useful as a 
theoretical idealization but is not relevant in practical cases, where the disorder 
distribution might be broad but remains normalizable. In that case, only the 
initial stage of damage accumulation resembles percolation, but eventually a 
catastrophic process is initiated and size effects, both in damage and strength 
are observed in agreement with experiments. 

The depinning approach is applicable instead to the weak disorder limit, 
single crack limit. Strong inhomogeneities would break up the single valued 
crack interface that is the starting assumption of the theory. According to the 
depinning scenario, the crack interface moves, eventually leading to the failure 
of the sample, when the stress intensity factor on the crack overcomes a critical 
threshold Kc- Since we deal here with a non-equilibrium phase transition, as 
in percolation, we do not expect to see size effects in Kc, besides the usual 
finite size corrections (i.e. > for L — > oo). Other possible size effects 
are hidden in the definition of the stress intensity factor. In other words, a 
depinning transition alone can not account for the observed size effects but 
can at best be used to describe the crack dynamics and morphology once crack 
growth has been initiated. 

7.2 Morphology of the fracture surface: roughness exponents 

We have discussed both the experimental (in Sec. l3.2|) . theoretical (in Sec. 15.31) 
and simulation-based (in Sec. 16. 4() crack roughness analysis. As we discussed 
extensively there are several possible definitions of the roughness exponent, de- 
pending on the experimental conditions (dimensionality of the sample, loading 
mode) and on the measured surface properties (in plane or out of plane, par- 
allel or perpendicular, local or global roughness). This variety of exponents 
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exponent 


exp. 


RFM 


DT-LR 


DT-SR 


C 2D local (i), 


0.6-0.7 


0.7 ± 0.1 






C 2D global (i) 




0.8 ± 0.1 






( 3D planar crack (ii) 


0.6 


fractal 


0.39 


0.63 


C± 3D local (iii) 


0.4-0.5 (ss) 0.8 (Is) 


0.4 


log 




C\\ 3D local (iii) 


0.6 


0.4 






C± 3D global (iii) 


1.2-1.3 


0.5 







Tabic 1. A summary of the values of the roughness exponent measured in experiments, in the RFM and 
in the depinning transition (DT) scenario, both with long-range (LR) and (SR) interactions. Refer to 
Sec. 12.6. ll for the definition of the exponents, to Sec. 13.21 for details about experiments, to Sec. Ib.^ for the 
RFM results and to Sec. l5.3l for the depinning transition. 



makes a direct comparison between theory and experiments quite involved. In 
order to simplify this task we report in Table ^ the values of the roughness 
exponents measured in experiments, as well as the predictions of the models. 

A direct comparison of the exponent values leads to two main conclusions: 
first, the exponents agree "roughly" in two dimensions (i.e. C — 0-7) - depend- 
ing on the case considered - but in three dimensions the numerical value is 
significantly smaller than the experimental one, Cioc — 0.4 and ^ ~ 0.5 in sim- 
ulations while Cioc — 0.8 and ^ ~ 1.2 in experiments. This discrepancy is quite 
striking given the large range over which ^joc can be observed in experiments. 
Second, The physical mechanism to crack roughening is still not understood 
even in the case of the RFM, not to mention other more realistic models and 
of course experiments. 

From the experimental point of view, presently it is not completely clear if 
the 3D fracture surface exponents are universal, and what kind of regimes can 
be found. The crucial lengthscales seem to be that of FPZ [89,92], ^fpz, and 
in some cases one related to voids formed ahead of the crack tip [95]: their 
density, and typical size. From the model point of view, it is crucial to underline 
here the main implication of the observed catastrophic failure occurring after 
peak-load and the lack of localization before this point: the roughness of the 
crack arises in the last abrupt "avalanche". This has important consequences 
for the modeling strategies, since quasi-static models such as the ones we are 
discussing here may become inappropriate when the crack grows unstable. It 
would be necessary to incorporate a more realistic dynamics, possibly including 
sound waves or other dynamical effects, or some process leading to residual 
strength. 

A further relevant theoretical point is the question of the universality of 
the roughness exponent, an hypothesis supported by ample experimental ev- 
idence. Early RFM simulations suggest that does not depend on the dis- 
order strength [349], while recent results indicate as well an independence of 
the lattice type [350]. It would be very interesting to consider in more de- 
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tail comparisons to/within truly vectorial (beam or bond-bending) models to 
see their roughness characteristics. This also holds for crack growth experi- 
ments from notches, in which the interaction of dominating crack with diffuse 
damage could be investigated. An important aspect stems from the connec- 
tion between fracture energy (or toughness) and the crack shape. It should be 
mentioned that the amplitude A of the roughness (as in oc AL^) usually 
increases as the prefactor of strength decreases if the disorder is varied. Thus 
we may conclude that "rough" does not necessarily imply "tough", to answer 
the original question posed 20 years ago in Ref. [88]. 



7.3 Crack dynamics: avalanches and acoustic emission 

The statistical properties of acoustic emission also remain to be explained, 
completely. Quasistatic lattice models qualitatively predict the correct behav- 
ior: a power law distribution of avalanche events, with an activity increasing 
as the failure point is approached. The scaling behavior is well captured by 
mean-field theory, as exemplified by fiber bundle models, with some deviations 
from the analytical predictions in the two-dimensional RFM (where r > 5/2), 
but almost perfect agreement in d = 3 and for the RSM. The apparent suc- 
cess of mean-field theory is probably related to the occurrence of a first-order 
transition in presence of long-range elastic interactions. In these conditions, 
the spinodal instability with the related mean-field scaling is accessible even 
in low dimensionality. From the point of view of the experimental comparison, 
however, lattice models typically predict a too large energy exponent /3. It is 
not known as to whether this problem results from any of the known possible 
caveats, such the brittle nature of the models, the dispersion of the acoustic 
waves and data in experiments, or the possibility that in reality one can not 
separate the time-scales of stress relaxation and loading, as assumed in the 
quasi-static models. For this reason we would like to call for more experiments 
in various loading conditions (tensile, creep, fatigue etc.), for varying materials 
so that the influence of plasticity, time-dependent response, and localization 
can be studied. 

Brittle fracture models are due to the first-order nature of the final fail- 
ure such that their predictability is inherently weak, though there have been 
(mainly experimental) attempts to show the opposite [121,151,361]. For the 
practically minded scientist, the main lessons are that brittle materials have 
a tendency to fail abruptly, exhibit critical-looking noise, and their fracture 
statistics depend crucially on the presence of large flaws and on the actual 
disorder. 
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7.4 From discrete models to damage mechanics 

Lattice models have been used to describe qualitatively the stress-strain curves 
of materials as noted earlier. Here, we have overviewed, mainly in the context 
of the RFM, the basic tenants of damage accumulation. The crucial questions 
are: how to average over samples, what kind of correlations arise in the damage, 
including localization, and in general what one learns from a representative 
volume element -idea. 

Large-scale simulations provide evidence for strong self-averaging for the D, 
damage variable, in the RFM. This is connected with the notion that while the 
damage accumulated becomes anisotropic, it still has only a finite correlation 
length ^ (both parallel and perpendicular to the load) up to the peak load. 
Meanwhile, as one takes the limit L ^ oo the strength decreases and actually 
the self- averaging gets stronger since L/^ remains large even up to the peak 
load. Again, it would be useful to test these conclusions in vectorial models. 
One question not addressed here in detail is the probability distribution of 
P{D, V), which is related to £,{V) via a simple application of the central limit 
theorem. Again, we would like to call for experiments, including looking at 
damage patterning and localization to also stimulate the further development 
of models. 



7.5 Concluding remarks and perspectives 

We think that there are substantial advances to be made within the confines of 
lattice fracture models with inelastic, or time-dependent effects incorporated. 
Even simulations of elastic-plastic fracture in this respect are still awaiting, 
though these have no rate-dependence. As was very briefiy mentioned in Sec. 
K14L there are statistical laws in fracture and deformation processes that are 
not simply "brittle" . For many of these cases we expect that models akin to 
the RFM and its generalizations will prove to be useful. 

There are some simple ingredients that "improved" models might contain 
and some simple developments they might lead to. For one, we would like to see 
the invention of lattice models that allow for true damage localization before 
the peak load. This would for instance produce damage dynamics comparable 
to the experimental "b- value" analysis of AE, to name but one issue. Similarly, 
one would need models that allow for stable crack growth beyond the peak 
load (F-controlled ensemble for equivalent RFM's). 

In all, theoretically it is also an interesting question how the "first-order" 
picture of fracture changes in the presence of other screening mechanisms for 
the long-range interactions. The implications of statistical fracture on other 
lengthscales, as in multi-scale modeling, are then dependent on what is the 
appropriate coarse-grained description. The resulting picture of the emerging 
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physics should be connected to what we have presented in this work. 
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Appendix A: Algorithms 

As mentioned earlier, fracture of disordered quasi-brittle materials is often 
studied using discrete lattice networks, wherein the damage is accumulated 
progressively by breaking one bond at a time until the lattice system falls 
apart. Each time a bond is broken, it is necessary to redistribute the stresses to 
determine the subsequent bond failure. Algebraically, this process of simulating 
fracture by breaking one bond at a time is equivalent to solving a new set of 
linear equations (Kirchhoff equations in the case of fuse models) 

A„x„ = b„, n = 0,1,2,..., (Al) 

every time a new lattice bond is broken. In Eq. ()A1|) . each matrix A„ is an 
N X N symmetric and positive definite matrix (the lattice conductance matrix 
in the case of fuse models and the lattice stiffness matrix in the case of spring 
and beam models), b„ is the A x 1 (given) applied nodal current or force 
vector, x„ is the A x 1 (unknown) nodal potential or displacement vector, 
and is the number of degrees of freedom (unknowns) in the lattice. The 
subscript n in Eq. ()A1|) indicates that A„ and b„ are evaluated after the n*'' 
bond is broken. The solution x„, obtained after the n*'' bond is broken, is used 
in determining the subsequent ((n -|- l)*'^) bond to be broken. 

Numerical simulations of fracture using large discrete lattice networks are 
often hampered due to the high computational cost associated with solving a 
new large set of linear equations every time a new lattice bond is broken. This 
becomes especially severe with increasing lattice system size, L, as the num- 
ber of broken bonds at failure, nj, follows a power-law distribution given by 
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Tif ~ 0{L^'^) in 2D and rif ~ 0(L^'^) in 3D. In addition, since the response of 
the lattice system corresponds to a specific reahzation of the random breaking 
thresholds, an ensemble averaging of numerical results over Neon fig configu- 
rations is necessary to obtain a realistic representation of the lattice system 
response. This further increases the computational time required to perform 
simulations on large lattice systems. 

A.l Fourier Acceleration Method 

Traditionally, iterative techniques based on preconditioned conjugate gradient 
(PCG) method have been used to simulate fracture using fuse networks (see 
Ref. [362] for an excellent review of iterative methods; see Ref. [363] for a re- 
view of multigrid method). However, large-scale numerical simulations using 
iterative solvers have often been hindered due to the critical slowing down asso- 
ciated with the iterative solvers as the lattice system approaches macroscopic 
fracture. That is, as the lattice system gets closer to macroscopic fracture, the 
condition number of the system of linear equations increases, thereby increas- 
ing the number of iterations required to attain a fixed accuracy. This becomes 
particularly significant for large lattices. As a remedy, Fourier accelerated PCG 
iterative solvers [351,364,365] have been suggested to alleviate the critical slow- 
ing down. The Fourier acceleration algorithm proposed in Refs. [364, 365] for 
solving the system of equations Ax = b can be expressed as shown in Algo- 
rithm 1. In this subsection, we do not explicitly write the subscript n to refer 
to the system of equations A„x„ = b„ obtained after the failure of n bonds. 
Instead, the subscript n is implicitly understood. 

Algorithm A.l 
1: Compute r^'^) = b — Ax^^) for some initial guess x^^) 
2: for /c = 1, 2, ... do 
3: Solve: Mz^'^-i) = r^^-^) 

4: xW = X^'^-^) + Z^''-!) 

6; Check convergence; continue if necessary 
7: end for 

In Algorithm 1, the matrix M is referred to as the preconditioner, and the 
superscript in the parenthesis refers to the iteration number. The steps in the 
Algorithm 1 can be combined into an iterative scheme as shown below 

^{k+i) ^ ^{k) ^ M-^j.{k) (A2) 

where r^'^^ = b — Ax*^'^'). Denoting the error e (or algebraic error) as e = x — x, 
where x is the exact solution such that Ax = b, the errors in the fc*^ and 
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{k + 1) iterations may be related as 

e(^'+i) = (I - M-iA)e('=) (A3) 

Consequently, after k number of relaxation sweeps, the error in the k^^ ap- 
proximation is given by e^*^) = Q'^e*^''), where Q = (I — M~^A). Choosing a 
particular vector norm and its associated matrix norm, it is possible to bound 
the error after k iterations by ||e'-'^)|| < ||Q||'^||e^'^^||. From this relation, it can 
be shown that the iteration associated with the matrix Q converges for all ini- 
tial guesses if and only if p{Q,) < 1, where p{Q) is the spectral radius defined 
as 

p(Q) =max,|Aj(Q)| (A4) 

and Aj(Q) denotes the j*^ eigenvalue of Q given by Aj(Q) = 1 — Aj(M^^A). 
Hence, convergence of Algorithm 1 is fast whenever the eigenvalues of M~^A 
are clustered around one, i.e., when is an approximation of the inverse 

of A. 

It should be noted that Eq. HA2|) is similar to the iterative algorithm based 
on the steepest descent (or alternatively, the modified Newton's method or the 
method of deflected gradients [366] ) , which is given by 

x{fc+i)^^(fc) + ^{fc)M-irW (A5) 

where a*^'^^ is a non-negative scalar chosen to minimize the standard quadratic 
function f{x^''^^^) defined as 

/(x) = ix^Ax - b^x (A6) 

The minimization of /(x) is equivalent to the solution of the system of equa- 
tions Ax = b. In this sense, a variety of iterative algorithms such as steepest 
descent and quasi-Newton methods [366] may be employed to solve Ax = b 
using the iterative scheme based on Eq. (|A5|) . For instance, = I results 

in a steepest descent algorithm, whereas if is the inverse of the Hessian 
of / then we obtain the Newton's method. In general, is chosen to be 

an approximation to the inverse of the Hessian. In addition, in order to guar- 
antee the iterative process defined by Eq. (|A5|) to be a descent method for 
small values of a^^\ it is in general required that M be positive definite [366]. 
For computational purposes, M is chosen such that its inverse can be com- 
puted cheaply so that the solution of Mz = r can be computed with least 
computational effort. 
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The Fourier acceleration algorithm proposed in Refs. [364, 365] chooses an 
ensemble averaged matrix A [367,368] as the preconditioner in Algorithm 1, 
where A{i,j) = r^'^™"'^^), r = \i — j\, the distance between the nodes i and 
j, and df and dw refer to the fractal dimension of the current-carrying back- 
bone and the random-walk dimension respectively. This ensemble averaged 
matrix A is a symmetric Toeplitz matrix and hence the preconditioned sys- 
tem Mz = r in Algorithm 1 can be solved in 0{N log N) operations by using 
FFTs of size N. Compared with the unconditioned CG methods, the Fourier 
accelerated technique based on ensemble averaged circulant preconditioner sig- 
nificantly reduced the computational time required to simulate fracture using 
random fuse networks [364, 365] . Although the Fourier accelerating technique 
significantly improves the condition number of the unconditioned system of 
equations, these methods still exhibit critical slowing down close to macro- 
scopic fracture. Consequently, the numerical simulations in the past were often 
limited to systems of sizes L = 128. 

Mathematically, the ensemble averaged matrix A is not the optimal circulant 
preconditioner in the sense that it does not minimize the norm ||I — C~^A||i7' 
over all non-singular circulant matrices C, and hence is expected to take more 
CG iterations than the optimal [369-372] and superoptimal [372, 373] circu- 
lant preconditioners. In the above description, || • Hi;' denotes the Frobenius 
norm [374]. It should be noted that ||I — C~^A||i7' determines the convergence 
rate of Algorithm 1 as evident in Eq. HA3|) . The discussion of optimal [369-372] 
and superoptimal [372, 373] circulant preconditioners is presented in detail in 
Section A. 4. In addition, the Fourier acceleration technique is not effective 
when fracture simulation is performed using central-force and bond-bending 
lattice models [365]. In the case of central-force and bond-bending lattice mod- 
els, each of the lattice sites has more than one degree of freedom. This results 
in a matrix Aq that is block-circulant. This explains why a naive circulant 
preconditioner is not effective in the case of central-force and bond-bending 
lattice models. A discussion on block-circulant preconditioners is presented in 
Section A. 4. 

A. 2 Naive Application of Sparse Direct Solvers 

Alternatively, since each of the matrices A„ for n = 0, 1, 2, . . . , is sparse, it is 
possible to use any of the existing sparse direct solvers [375,376] to solve each 
of the system of equations formed by A„. It should be noted that the Cholesky 
factorization of A„ for any particular n can be done efficiently. However, fac- 
torization of each of the A„ matrices for ~ 0{L^'^) (or ~ 0{Lp' '^) in 
3D) number of times is still a computationally daunting task. For example, 
in the case of a 2D triangular lattice system of size L = 128, the algorithm 
based on factorizing each of the A„ matrices using supernodal sparse direct 
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solver [377] took on an average 8300 seconds compared to an average of 7473 
seconds taken by the optimal circulant preconditioner based CG algorithm. It 
should however be noted that, in general, the optimal circulant preconditioner 
based CG algorithm used in the above comparison is much faster than the 
Fourier accelerated CG algorithm described in Ref. [365]. However, as the fol- 
lowing description elaborates, the computational advantages of sparse direct 
solvers are not fully realized when the algorithm is based on factorizing each 
of the An matrices. 

A. 3 Multiple-rank sparse Cholesky downdating algorithm 

An important feature of fracture simulations using discrete lattice models is 
that, for each n = 0, 1, 2, . . ., the new matrix A„_|_i of the lattice system after 
the (n + l)*'^ broken bond is equivalent to a rank-p downdate of the matrix 
An [187,378]. Mathematically, in the case of the fuse and spring models, the 
breaking of a bond is equivalent to a rank-one downdate of the matrix A„, 
whereas in the case of beam models, it is equivalent to multiple-rank (rank-3 
for 2D, and rank-6 for 3D) downdate of the matrix A„. Thus, an updating 
scheme of some kind is therefore likely to be more efficient than solving the new 
set of equations formed by Eq. ()A1|) for each n. One possibility is to use the 
Shermon-Morrison- Woodbury formula [374] to update the solution vector x^+i 
directly knowing the sparse Cholesky factor of A^, for any m < n. That 
is, given the Cholesky factorization of A^ for m = 0, 1, 2, • • •, it is possible 
to obtain a direct updating of the solution x„+i for n = m,m + l,m + 2, ■ ■ ■ 
based on p = (n — m) saxpy vector updates [187,378]. The resulting algorithm 
is equivalent to a dense matrix update algorithm proposed in the Appendix C 
of Ref. [187]. In Ref. [187], it has been shown that for two-dimensional fracture 
simulations, any of these sparse solver updating algorithms are significantly 
superior to the traditional Fourier accelerated CG iterative schemes. The best 
performing sparse solver for two-dimensional fracture simulations using the 
discrete lattice models (fuse, spring and beam models) is based on the multiple- 
rank sparse Cholesky downdating algorithm, which is briefly described in the 
following. 

Consider the Cholesky factorizations 

PA„P* = l^nK (A7) 

for each n = 0, 1, 2, • • •, where P is a permutation matrix chosen to preserve 
the sparsity of L„. Since the breaking of bonds is equivalent to removing the 
edges in the underlying graph structure of the matrix A„, for each n, the 
sparsity pattern of the Cholesky factorization L„+i of the matrix A„_|_i must 
be a subset of the sparsity pattern of the Cholesky factorization L„ of the 
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matrix A„. Hence, for all n, the sparsity pattern of L„ is contained in that of 
Lq. That is, denoting the sparsity pattern of L by C, we have 

jO-m^^n ^ m<n (A8) 

Therefore, we can use the sparse Cholesky factorization downdate algorithm 
of Davis and Hagcr [379, 380] to successively downdate the Cholesky factor- 
izations L„ of A„ to L„+i of A„+i, i.e., L„ L„_|_i for n = 0, 1, 2, • • •. Since 
i^n 2 -Cn+i, it is necessary to modify only a part of the non-zero entries of 
L„ in order to obtain L„ ^n+i- This results in a significant reduction in 
the computational time. In addition, the rank-p downdating vectors are very 
sparse with only a few non-zero entries (at most two for fuse models; six for 
spring models; and twelve for beam models). Consequently, the computational 
effort involved in downdating the Cholesky factors using the sparse downdat- 
ing vectors is once again significantly reduced [187,378]. 

A. 4 Optimal and superoptimal circulant preconditioners 

Although the algorithms based on sparse direct solvers with multiple-rank 
Cholesky downdating schemes achieve superior performance over iterative 
solvers in 2D lattice simulations (see Section A. 7), the memory demands 
brought about by the amount of fill-in during sparse Cholesky factorization 
poses a severe constraint over the usage of sparse direct solvers for 3D lat- 
tice simulations. Hence, iterative solvers are in common use for large-scale 3D 
lattice simulations. The performance of iterative solvers depends crucially on 
the condition number and clustering of the eigenvalues of the preconditioned 
system. In general, the more clustered the eigenvalues are, the faster the con- 
vergence rate is. Hence, the main focus of any iterative solution technique is 
choosing an optimum preconditioncr. 

The main observation behind developing preconditioners for the iterative 
schemes is that the operators on discrete lattice network result in a circulant 
block structure. For example, the Laplacian operator (Kirchhoff equations in 
the case of random fuse model) on the initial uniform grid results in a Toeplitz 
matrix Aq. Hence, a fast Poisson type solver with a circulant preconditioncr 
can be used to obtain the solution in 0{N logN) operations using FFTs of 
size N, where N denotes the number of degrees of freedom. However, as the 
lattice bonds are broken successively, the initial uniform lattice grid becomes a 
diluted network. Consequently, although the matrix Aq is Toeplitz (also block 
Toeplitz with Toeplitz blocks) initially, the subsequent matrices A„, for each 
n, are not Toeplitz matrices. However, depending on the pattern of broken 
bonds, A„ may still possess block structure with many of the blocks being 
Toeplitz blocks. 
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The optimal circulant preconditioner c(A) [369-372] of a matrix A is defined 
as tfie minimizer of ||C — A||i? over all x circulant matrices C. Given a 
matrix A, the optimal circulant preconditioner c(A) is uniquely determined 
by 

c(A) =F*(5(FAF*)F (A9) 

where F denotes the discrete Fourier matrix, 5 (A) denotes the diagonal matrix 
whose diagonal is equal to the diagonal of the matrix A, and * denotes the 
adjoint (i.e. conjugate transpose). It should be noted that the diagonals of 
FAF* represent the eigenvalues of the matrix c(A) and can be obtained in 
0{N log N) operations by taking the FFT of the first column of c(A). The 
first column vector of T. Chan's optimal circulant preconditioner matrix that 
minimizes the norm ||C — A||i? is given by 

1 ^ 

^ 77 X] ""i.O-i+i) mod N (AlO) 
i=i 

If the matrix A is Hermitian, the eigenvalues of c(A) are bounded below and 
above by 

(C(A)) < XmaxidA)) < A 

max 

(A) (All) 

where Amm(") and Xmaxi') denote the minimum and maximum eigenvalues, 
respectively. Based on the above result, if A is positive definite, then the 
circulant preconditioner c(A) is also positive definite. The computational cost 
associated with the solution of the preconditioned system c(A)z = r is the 
initialization cost of nnz(A) for setting the first column of c(A) using Eq. 
(|A10|) during the first iteration, and 0{N log N) during every iteration step. 

The superoptimal circulant preconditioner t{A) [373] is based on the idea of 
minimizing the norm ||I — C~^A||i;' over all nonsingular circulant matrices C. 
In the above description, t{A) is superoptimal in the sense that it minimizes 
||I — C^^A||i?, and is equal to 

t(A) = c(AA*)c(A)-^ (A12) 

The preconditioner obtained by Eq. (|A12|) is also positive definite if A itself is 
positive definite. Although the preconditioner t(A) is obtained by minimizing 
the norm ||I — C~^A||i;', the asymptotic convergence of the preconditioned 
system is same as c(A) for large A^ system. 
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A. 5 Block- circulant preconditioner 

Alternatively, block-circulant preconditioners may also be used as precondi- 
tioners since many of the block matrices in for n > may still retain the 
circulant or Toeplitz property. Let the matrix A„ be partitioned into r-by-r 
blocks such that each block is an s-by-s matrix. That is, = rs. Since each of 
the blocks of A„ are Toeplitz (or even circulant), the block-circulant precon- 
ditioner is obtained by using circulant approximations for each of the blocks 
of An- It is the minimizcr of ||C — A„||i,^ over all matrices C that are r-by-r 
block matrices with s-by-s circulant blocks. In addition, we have 

max max (An) (A13) 

In particular, if A„ is positive definite, then the block-preconditioner cb(A„) 
is also positive definite. In general, the average computational cost of using 
the block-circulant preconditioner per iteration is 0{rs log s) + delops, where 
delops represents the operational cost associated with solving a block-diagonal 
matrix with r x r dense blocks. For 2D and 3D discrete lattice network with 
periodic boundary conditions in the horizontal direction, this operational cost 
reduces significantly. The reader is referred to [187, 381] for further details on 
optimal and block-circulant preconditioners for discrete lattice networks. 

A. 6 Summary of algorithms 

To summarize, we use the following two algorithms based on sparse direct 
solvers for modeling fracture using discrete lattice systems. 

• Solver Type A: Given the factorization L„j of A„j, the rank-1 sparse 
Cholesky modification is used to update the factorization L„+i for all sub- 
sequent values of n = m, m + 1, • • •. Once the factorization L„_|_i of A^+i is 
obtained, the solution vector x„_|_i is obtained from L„+iL^^^x„+i = b^+i 
by two triangular solves [187,378]. 

• Solver Type B: Assuming that the factor L^, of matrix A^ is available after 
the m*'* fuse is burnt, the solution x„_|_i after the (n -|- 1)*'* fuse is burnt, for 
n = m, m + 1, ■ ■ ■ , m + maxupd, is obtained through p = [n — m) saxpy vector 
updates using an algorithm based on Shermon-Morrison- Woodbury [374] 
formula. Either a completely new factorization of the matrix A is performed 
or a multiple-rank downdate of the factorization — ^ "^m+maxupd+i is 
performed every maxupd steps [187,378]. 

In addition to the above two algorithms, the simulations are carried out 
using the iterative solvers based on optimal circulant and block-circulant pre- 
conditioners [187,381]. 
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A. 7 Performance comparison 

For two-dimensional lattice systems, we choose a triangular lattice topology, 
whereas a cubic lattice topology is used in three-dimensional cases. The ran- 
dom thresholds fuse model (rank-one downdate with = 1) is used as a basis 
of performance comparison to evaluate the numerical efficiency of various algo- 
rithms. The algorithmic framework for central-force (spring) and beam models 
is presented in Refs. [187,378]. 

In the numerical simulations using solver types A and B based on sparse 
direct solvers, the maximum number of vector updates, maxupd, is chosen 
to be a constant for a given lattice size L. We choose maxupd = 25 for 
L = {4, 8, 16, 24, 32}, maxupd = 50 for L = 64, and maxupd = 100 for 
L = {128,256,512}. For L = 512, maxupd is limited to 100 due to mem- 
ory constraints. By keeping the maxupd value constant, it is possible to com- 
pare realistically the computational cost associated with different solver types. 
Moreover, the relative CPU times taken by these algorithms remain the same 
even when the simulations are performed on different platforms. 

The convergence rate of an iterative solver depends on the quality of starting 
vector. When PCG iterative solvers are used in updating the solution x„ 
x„+i, we use either the solution from the previous converged state, i.e., x„, or 
the zero vector, 0, as the starting vectors for the CG iteration of x„+i. The 
choice of using either x„ or as the starting vectors depends on which vector 
is closer to hn+i in 2-norm. That is, if ||b„+i — A„x„||2 < ||b„+i||2, then x„ is 
used, else is used as the starting vector for the x„+i CG iteration. A relative 
residual tolerance e = 10~^^ is used in CG iteration. 

The numerical simulations are carried out on a single 1.3 GHz IBM Power4 
processor. Tables Al and A2 present the CPU times taken for one configuration 
(simulation) using the sparse direct solver algorithms A and B, respectively. 
These tables also indicate the number of configurations. Neon fig, over which 
ensemble averaging of the numerical results is performed for performance eval- 
uation. The CPU times taken by the iterative solvers are presented in Tables 
A3-A4. For iterative solvers, the number of iterations presented in Tables A3- 
A4 denote the average number of total iterations taken to break one intact 
lattice configuration until it falls apart. In the case of iterative solvers, some 
of the simulations for larger lattice systems were not performed either because 
they were expected to take larger CPU times or the numerical results do not 
influence the conclusions drawn in this study. The results presented in Tables 
A1-A4 indicate that for 2D lattice systems, the multiple-rank Cholcsky down- 
date based sparse direct solver algorithms are clearly superior to the optimal 
and block-circulant preconditioned CG iterative solvers. In particular, for the 
lattice size L = 128, the solver type A took 212.2 seconds, which should be 
compared with the 7473 and 8300 seconds taken respectively by the optimal 
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Table Al. Solver Type A (2D Triangular Lattice) 



Size 


CPU(sec) 


Simulations 


32 


0.592 


20000 


64 


10.72 


4000 


128 


212.2 


800 


256 


5647 


96 


512 


93779 


16 



Table A2. Solver Type B (2D Triangular Lattice) 



Size 


CPU(sec) 


Simulations 


32 


0.543 


20000 


64 


11.15 


4000 


128 


211.5 


800 


256 


6413 


96 



Table A3. Block Circulant PCG (2D Triangular Lattice) 



Size 


CPU(sec) 


Iter 


Simulations 


32 


10.00 


11597 


20000 


64 


135.9 


41207 


1600 


128 


2818 


147510 


192 


256 


94717 




32 



Table A4. Optimal Circulant PCG (2D Triangular Lattice) 



Size 


CPU(sec) 


Iter 


Simulations 


32 


11.66 


25469 


20000 


64 


173.6 


120570 


1600 


128 


7473 


622140 


128 



circulant preconditioner and the naive sparse direct solution technique based 
on factorizing each of the A„ matrices. 



For large 3D lattice systems, the advantages exhibited by the multiple-rank 
Cholesky downdating algorithms in 2D simulations vanish due to the amount 
of fill-in during Cholesky factorization. Tables A5-A8 present the CPU times 
taken for simulating one-configuration using the block circulant, optimal cir- 
culant, un-preconditioned, and the incomplete Cholesky iterative solvers, re- 
spectively. It should be noted that for large 3D lattice systems (e.g., L = 32), 
the performance of incomplete Cholesky preconditioner (see Table A8) is sim- 
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Table A5. Block Circulant PCG (3D Cubic Lattice) 



Size 


CPU(sec) 


Iter 


Simulations 


10 


16.54 


16168 


40000 


16 


304.6 


58756 


1920 


24 


2154 


180204 


256 


32 


12716 


403459 


128 


48 


130522 


1253331 


32 


64 


1180230 




11 



Table A6. Optimal Circulant PCG (3D Cubic Lattice) 



Size 


CPU(sec) 


Iter 


Simulations 


10 


15.71 


27799 


40000 


16 


386.6 


121431 


1920 


24 


2488 


446831 


256 


32 


20127 


1142861 


32 


48 


233887 


4335720 


32 



Table A7. Un-preconditioncd CG (3D Cubic Lattice) 



Size 


CPU(sec) 


Iter 


Simulations 


10 


5.962 


48417 


40000 


16 


119.4 


246072 


3840 


24 


1923 


1030158 


256 


32 


16008 


2868193 


64 



Table A8. Incomplete Cholosky PCG (3D Cubic Lattice) 



Size 


CPU(scc) 


Iter 


Simulations 


10 


5.027 


8236 


40000 


16 


118.1 


42517 


3840 


24 


1659 


152800 


512 


32 


12091 


122113 


(34 



ilar to that of block-circulant preconditioner (see Table A5), even though the 
performance of incomplete Cholcsky preconditioner is far more superior in the 
case of 2D lattice simulations [187]. Indeed, for 3D lattice systems of sizes 
L = 48 and L = 64, the CPU time taken by the sparse direct solver types A 
and B is much higher than that of block-circulant PCG solver. 

To summarize, the numerical simulation results presented in Tables A1-A4 
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for 2D triangular random fuse networks indicate that sparse direct solvers 
based on multiple-rank sparse Cholesky downdating are superior to the com- 
peting iterative schemes using block-circulant and optimal circulant precon- 
ditioners. On the other hand, for 3D lattice systems, the block-circulant prc- 
conditioner based CG solver exhibits superior performance (for system sizes 
L > 32) over the sparse direct solvers and the related incomplete Cholesky 
preconditioned CG solvers. In addition, for both 2D and 3D lattice systems, 
the block circulant preconditioned CG is superior to the optimal circulant pre- 
conditioned PCG solver, which in turn is superior to the Fourier accelerated 
PCG solvers used in Refs. [364,365]. 
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